83

OK, there is good news and there is bad news. In the current version 10 there is no way to do this directly. That's the bad news. The good news is that finite element framework used within NDSolve is exposed and documented; for maximum "hackability" convenience. Let's start with a region that @MarkMcClure would consider interesting. We load our favorite ...


66

The only reason I am attempting to answer this is to perhaps get a Reversal badge. There you go... We will go slowly and this answer is the basis for what comes next. Let's start with two dimensions. You'll see why. We create a rectangular region: Needs["NDSolve`FEM`"] mesh = ToElementMesh[FullRegion[2], {{0, 5}, {0, 1}}, "MeshOrder" -> 1, "...


60

I've encapsulated the code of the mysterious user21 into a helmholzSolve command. The code is at the end of this post. It adds very little to user21's code but it does allow us to examine multiple examples quite easily, though it has certainly not been tested extensively and could be improved quite a lot I'm sure. It should be called as follows: {ev,if,...


54

Background Details about multigrid solvers can be found in this pretty neat script by Volker John. That's basically the source from which I drew the information to implement the V-cycle solver below. In a nutshell, a multigrid solver builds on two ingredients: A hierarchy of linear systems (with so-called prolongation operators mapping between them) and a ...


42

I can give a specific opinion on the comparison. The specificity is that I am not at all a specialist in numerical methods, but rather a user in this area, a physicist who needs to solve equations. First about strong features of Comsol with respect to Mathematica (MMA). The strongest feature of Comsol is that it has a nonlinear solver, while MMA still ...


40

OK, let me come straight - I did not read your question much beyond the title and this post will not address the specific issues you raise in your question. As an alternative I'll show how to use the low level FEM functionality to code up a non-linear Navier-Stokes solver. The documentation explains the details about the low level FEM programming ...


37

Version 11 has both symbolic and numeric eigensolvers, see here for an overview Here is a slightly different way to do it. We write a function that converts any PDE (1D/2D/3D) into discretized system matices: Needs["NDSolve`FEM`"] PDEtoMatrix[{pde_, Γ___}, u_, r__, o : OptionsPattern[NDSolve`ProcessEquations]] := Module[{ndstate, feData, sd, bcData, ...


36

After some self-learning of finite element method (FEM), I'm able to write a solver for OP's question without NDSolve now. The following code is highly inspired by ruebenko's excellent answer here. To understand the code, you must have some basic understanding for FEM, which anyone with a conscience won't say it's trivial 囧. The best introduction for FEM I ...


26

Introduction I think there are several questions on this site about ODEs of the form $$(x-a)^2 u''(x) = F(x,u,u')$$ with an initial condition at $x=a$. There is no general guarantee that solutions exist over an interval $(a,b]$, but sometimes it is possible as in this case. Outline We transform the equation $u''(x) = F(x,u,u')$ over the infinite interval ...


24

Here are a couple of suggestions. First let's look at the mesh generation with is documented with ToElementMesh and in the mesh generation tutorial. Needs["NDSolve`FEM`"]; x2 = 2; y2 = 1; r = 1/8; reg = ImplicitRegion[ 0 <= x <= x2 && 0 <= y <= y2 && (x - x2/2)^2 + (y - y2/2)^2 >= r^2, {x, y}]; For the mesh generation ...


24

OK, there are a few things going on here. Let me explain them in turn. First, as the message suggests, this should be written in Inactive form (we'll get to the why later). If you click on the three dots in front of the error message and follow the link to the reference page you will find some information on this error message. To write the equation in ...


22

As of v10.3, this is now possible using the built-in commands DEigensystem and NDEigensystem functions. These work pretty much as you hope they would, along the lines of DEigensystem[ℒ[u[x, y, …]], u, {x, y, …} ∈ Ω, n] gives the n smallest magnitude eigenvalues and eigenfunctions for the linear differential operator ℒ over the region Ω. Hat-tip to ...


22

I have not specifically used COMSOL, but can speak to other FEA and CFD solvers that I use including ANSYS, LS_DYNA, ABAQUS, Fluent, CFX, etc. If you are doing engineering simulation of anything other than very simple problems (like the kinds you program in college classes) then a commercial package is vital. Speed is probably the least important ...


21

Lets start by making a symbolic representation of maze, with BooleanRegion using graphic primitives. That will be our basis for discretization by any method. Number of vertices of graph should be small to make experimentation easier. (* "IGraphM" package downloaded from https://github.com/szhorvat/IGraphM *) Get["IGraphM`"] n = 3; (* Number of vertices per ...


20

first part..i had lying around.. poly = Random[Real, {1, 2}] {Cos[#], Sin[#]} & /@ Sort[Table[Random[Real, {0, 2 Pi}], {5}]] isLeft[P2_, {P0_, P1_}] := -Sign@Det@{P2 - P0, P1 - P0}; pinpoly[p_, poly_] := Module[{ed},(*winding rule*) ed = Partition[Append[poly, poly[[1]]], {2}, 1]; Count[ed,pr_ /; (pr[[1, 2]] <= p[[2]] < pr[[2, 2]] &&...


20

Based on the tutorial here QuadElement meshes behave exactly the same as TriangleElement meshes, with the exception that, for linear quad elements, four incidents per element are needed, and, for quadratic elements, eight incidents per element are needed. here is another possible answer: Needs["NDSolve`FEM`"] rectangleA = {{0, 4}, {8, 8}}; rectangleB = ...


19

Using V10's new FEM functionality, this problem can be solved as follows << NDSolve`FEM`; omega = ImplicitRegion[x^6 + y^4 <= 1, {x, y}]; mesh = ToElementMesh[omega, "MaxCellMeasure" -> {"Area" -> 0.005, "Length" -> 0.1}]; gamma = DirichletCondition[u[t, x, y] == 0, x^6 + y^4 == 1]; u = NDSolveValue[ {Derivative[1, 0, 0][u][t, x, y] ==...


19

A useful feature that I regularly use in COMSOL and would like to be able to use in Mma is the "AdaptiveMeshRefinement" (as it is called in COMSOL). This means that COMSOL makes a mesh. With this mesh, it solves the problem. Then it evaluates a function that characterizes the steepness of the solution. Typically, it is the gradient of the solution squared,...


18

Table[drawtriangulation[mesh @@ example, First@example, AspectRatio -> Automatic], {example, {circle, circle34, ellipseeye}}] // GraphicsRow Calculating specifications for these examples: (* distance function, bounding box, fixed points, number of initial points, max iterations, min triangle quality *) circle = {Sqrt[#1^2 + #2^2] - 1. &, {{-...


18

EDIT from the author of this answer, 20 March 2016 This old answer was preceding Mathematica Version 10 which affords the finite Element Method. So this answer is now obsolete. See other answers below. The only interesting information that is remaining is that the Method of lines associated with the TensorProductGrid discretization can eventually handle ...


18

This slightly modified function Needs["NDSolve`FEM`"]; helmholzSolve3D[g_, numEigenToCompute_Integer, opts : OptionsPattern[]] := Module[{u, x, y, z, t, pde, dirichletCondition, mesh, boundaryMesh, nr, state, femdata, initBCs, methodData, initCoeffs, vd, sd, discretePDE, discreteBCs, load, stiffness, damping, pos, nDiri, ...


18

I have used both Mathematica and Comsol Multiphysics (formerly FEMLAB) for decades to solve practical problems. Mathematica is my go-to mathematics application because of versatility. I go to Comsol when: 1. There is a geometry more complicated than a box, cylinder, or sphere. Comsol uses consecutive solid geometry to build up shapes including the ...


17

Too long for a comment: OK, first set up your system that you only have non periodic BCs. Then look at the Finite Element programming tutorial and use NDSolve ProcessEquations and follow the steps until the call to DiscretizePDE and DiscretizeBoundaryConditions. At this point you can extract the system matrices. Deploy the (non periodic) boundary conditions....


17

In Version 12.0 this has become much easier, as now there is a nonlinear finite element solver. Here is the code: Set up the region. mu = 10^-3; rho = 1; l = 2.2; h = 0.41; region = RegionDifference[Rectangle[{0, 0}, {l, h}], Disk[{1/5, 1/5}, 1/20]]; Setup the equation. op = { rho*{u[x, y], v[x, y]}.Inactive[Grad][u[x, y], {x, y}] + Inactive[...


17

There are several part to you question and I'll try to tackle them one by one (not sure I'll be able to do this in one go) How to calculate stresses Here is the model problem. L = 1; h = 0.125; (*Shear stress on beam*) ss = 5; reg = Rectangle[{0, -h}, {L, h}]; mesh = ToElementMesh[reg]; materialParameters = {Y -> 10^3, \[Nu] -> 33/100}; We use a ...


17

Here is a different approach. The idea is to first generate a second order triangle mesh. Next a center coordinate is added in every triangle. Then we split every triangle into three first order quads making use of the newly added center coordinate and the initial second order triangle incidents. So we add a node in the center of the element and create the ...


17

3D Example The problem with direct solvers is that starting in 3 dimensions, their performance for dealing with matrices stemming from PDEs drops rapidly. This is why I wanted to show at least one 3-dimensional example. As there is no immediate analogon for Loop subdivision of tetrahedral meshes, I use hexahedral meshes instead. Preparation Here are some ...


17

Here's an approach that unions the primitives as rasters, meshes, and smooths. Data from OP: showmaze = Uncompress[FromCharacterCode @@ ImageData[Import["https://i.stack.imgur.com/XVJcP.png"], "Byte"]]; Primitives with an extended inlet and outlet: prims = CapsuleShape @@@ Cases[showmaze, _Cylinder, ∞]; prims = prims /. {{5., 5., 5.} -> {5.5, 5., 5.}, ...


Only top voted, non community-wiki answers of a minimum length are eligible