4
$\begingroup$

I have only recently been introduced to Mathematica's(v10.0) FEM capabilities. I understand that for solving PDEs on non-uniform shapes via NDSolve, Mathematica uses FEM. I have been able to extract the mesh from some examples shown and I notice that the mesh used by Mathematica is always ordered with equal element sizes. Is there some way I could force NDSolve to choose a combination of quad and tri elements instead of ONLY quad elements or tri elements to solve a PDE? Also is mesh refinement possible (dense mesh in some areas and coarse mesh in others)?

Here's a MWE that I would like to use: It is a notebook file. Forgive the link instead of pasting input text here: the input text was quite garbled and didn't copy back into a Mathematica.


update (to copy code from notebook here)

Ω = 
  ImplicitRegion[ ! (x - 5)^2 + (y - 5)^2 <= 3^2, {{x, 0, 5}, {y, 0, 
     10}}]; 
RegionPlot[Ω, AspectRatio -> Automatic]

op = -Laplacian[u[x, y], {x, y}] - 20; 

Subscript[Γ, 
   D] = {DirichletCondition[u[x, y] == 0, x == 0 && 8 <= y <= 10], 
   DirichletCondition[u[x, y] == 100, 
         (x - 5)^2 + (y - 5)^2 == 3^2]}; 

uif = NDSolveValue[{op == 0, Subscript[Γ, D]}, u, 
  Element[{x, y}, Ω], 
     Method -> {"PDEDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.3}, 
      "IntegrationOrder" -> 5}}]

Quiet[ContourPlot[uif[x, y], Element[{x, y}, Ω], 
  ColorFunction -> "Temperature", AspectRatio -> Automatic]]
Plot3D[uif[x, y], Element[{x, y}, Ω], 
 ColorFunction -> "Temperature", AspectRatio -> Automatic, 
 Mesh -> False]

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];

{mesh} = InterpolatingFunctionCoordinates[uif]
mesh["Wireframe"]

Mathematica graphics

$\endgroup$
  • $\begingroup$ fyi, I copied the code from the notebook. To do this in the future, convert the cells to InputForm and then copy as text and paste here as code. $\endgroup$ – Nasser May 18 '15 at 0:24
  • 2
    $\begingroup$ There are some examples of mixed meshes and difference densities in the documentation for ToElementMesh. See "Scope" and the option "MeshRefinementFunction". $\endgroup$ – Michael E2 May 18 '15 at 3:06
  • 2
    $\begingroup$ Could you give a hint, of where do you need the quad and where the triangle mesh? Do you need the mesh refined everywhere, of you want it to be refined locally? If locally, where? $\endgroup$ – Alexei Boulbitch May 18 '15 at 8:19
  • 2
    $\begingroup$ No, I think you can do it. However, it doubt it will happen automatically if you use ToElement mesh on a region or a BoundaryMesh. You may have to do it by hand, at least in part, as it is shown in the docs. If I get a chance, I'll try to experiment. $\endgroup$ – Michael E2 May 18 '15 at 18:24
  • 1
    $\begingroup$ @drN, you can have both tri and quad elements in a 2D mesh. You'd either create this manualy or import a mesh from another source; There is currently (V10.1) no functionality to generate pure quad (or mixed quad and tri) element meshes for non rectangular regions automatically. (I think that would be an interresting project, to write some Mathematica code that uses a triangle mesh and onverts it to a quad dominant mesh). For the refinement, as others pointed out the documentation has examples of how to write a refineQ function. $\endgroup$ – user21 Jun 5 '15 at 12:30
4
$\begingroup$

It's a bit late but here is a way to do it. The idea is to generate a quad mesh from a part of the boundary. The interface to the triangle mesh will be used from the quad mesh and we will generate the triangle mesh without adding Steiner points. This will help us to merge the meshes later into a single mesh.

r1 = Rectangle[{0, 0}, {3/2, 10}];
r2 = RegionDifference[Rectangle[{3/2, 0}, {5, 10}], Disk[{5, 5}, 3]];
Needs["NDSolve`FEM`"]
(leftMesh = ToElementMesh[r1, "MeshOrder" -> 1])["Wireframe"]

enter image description here

Next, we generate a boundary mesh of r2

rightBoundaryMesh = ToBoundaryMesh[r2, "MaxBoundaryCellMeasure" -> 0.1];

The idea is to select the right most edge from the quad mesh and the rest form the right boundary mesh to form a new boundary mesh with the spacing of the boundary elements of the quad mesh on the left edge of the right boundary.

coordinates = 
  Join[Select[leftMesh["Coordinates"], (#[[1]] == 3/2) &], 
   Select[rightBoundaryMesh["Coordinates"], (#[[1]] > 3/2) &]];
incidents = Partition[FindShortestTour[coordinates][[2]], 2, 1];
rightBoundaryMeshNew = ToBoundaryMesh["Coordinates" -> coordinates, 
   "BoundaryElements" -> {LineElement[incidents]}];

Now, we generate a boundary mesh of the right part that uses the boundary element spacing of the quad mesh for the left part. The key point is to use "SteinerPoints" -> False such that no new nodes are introduced on the edges.

rightMesh = 
  ToElementMesh[rightBoundaryMeshNew, "MeshOrder" -> 1, 
   "SteinerPoints" -> False];

rightMesh["Wireframe"]

enter image description here

Note, how on the left edge the quad boundary element spacing is used while on the three remaining edges a "MaxBoundaryCellMeasure" -> 0.1 is used (for illustration).

The last step is to glue the meshes together:

leftCoordinates = leftMesh["Coordinates"];
rightCoordinates = rightMesh["Coordinates"];
mesh = ToElementMesh[
   "Coordinates" -> Join[rightCoordinates, leftCoordinates], 
   "MeshElements" -> 
    Flatten[{rightMesh["MeshElements"], 
      QuadElement /@ (ElementIncidents[leftMesh["MeshElements"]] + 
         Length[rightCoordinates])}]];
mesh["Wireframe"]

enter image description here

This is a first order mesh. You could generate a second order mesh with

MeshOrderAlteration[mesh, 2]

But then you'd need to adjust the mid side nodes of the curved part. Which is certainly doable.

Solve and visualize a PDE:

op = -Laplacian[u[x, y], {x, y}] - 20;
Subscript[\[CapitalGamma], 
   D] = {DirichletCondition[u[x, y] == 0, x == 0 && 8 <= y <= 10], 
   DirichletCondition[u[x, y] == 100, (x - 5)^2 + (y - 5)^2 == 3^2]};
uif = NDSolveValue[{op == 0, Subscript[\[CapitalGamma], D]}, u, 
   Element[{x, y}, mesh]];
Show[
 Plot3D[uif[x, y], Element[{x, y}, mesh], 
  ColorFunction -> "Temperature", AspectRatio -> Automatic, 
  Mesh -> False, Boxed -> False, Axes -> None],
 Graphics3D[{Directive[{EdgeForm[Gray], FaceForm[]}], 
   ElementMeshToGraphicsComplex[mesh, 
    "CoordinateConversion" -> (Join[#, ConstantArray[{0.}, Length[#]],
         2] &)]}]
 ]

enter image description here

$\endgroup$
  • $\begingroup$ +1, I didn't know about the usefulness of "SteinerPoints" -> False before. $\endgroup$ – Pinti Apr 24 at 9:24
  • $\begingroup$ @Pinti, the usage shown here is, as far as I know, the primary usage: generate meshes that do not modify the input boundary. $\endgroup$ – user21 Apr 24 at 9:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.