3
$\begingroup$

I would like to use ElementMarker in a DirichletCondition on a structured quad mesh, but I am receiving the error:

NDSolveValue::bcnop: No places were found on the boundary where ElementMarker==2 was True, so DirichletCondition[u==1,ElementMarker==2] will effectively be ignored.

Here is an example of a small mesh where ElementMarker==2 appears to be on the bottom boundary.

Needs["NDSolve`FEM`"]
bounds = <|"inlet" -> 1, "hot" -> 2, "outlet" -> 3|>;
regs = <|"solid" -> 10, "fluid" -> 20, "interface" -> 15|>;
crd = {{0.`, 0.`}, {0.`, 0.4002986944615309`}, {0.`, 
    0.8326451978805829`}, {0.`, 1.2996052494743657`}, {0.`, 
    1.803950000871885`}, {0.`, 2.348672461377994`}, {0.`, 
    2.9370052598409973`}, {0.`, 3.5724398285307286`}, {0.`, 
    4.258747122872904`}, {0.`, 5.`}, {1.`, 0.`}, {1.`, 
    0.4002986944615309`}, {1.`, 0.8326451978805829`}, {1.`, 
    1.2996052494743657`}, {1.`, 1.803950000871885`}, {1.`, 
    2.348672461377994`}, {1.`, 2.9370052598409973`}, {1.`, 
    3.5724398285307286`}, {1.`, 4.258747122872904`}, {1.`, 5.`}, {2.`,
     0.`}, {2.`, 0.4002986944615309`}, {2.`, 
    0.8326451978805829`}, {2.`, 1.2996052494743657`}, {2.`, 
    1.803950000871885`}, {2.`, 2.348672461377994`}, {2.`, 
    2.9370052598409973`}, {2.`, 3.5724398285307286`}, {2.`, 
    4.258747122872904`}, {2.`, 5.`}, {3.`, 0.`}, {3.`, 
    0.4002986944615309`}, {3.`, 0.8326451978805829`}, {3.`, 
    1.2996052494743657`}, {3.`, 1.803950000871885`}, {3.`, 
    2.348672461377994`}, {3.`, 2.9370052598409973`}, {3.`, 
    3.5724398285307286`}, {3.`, 4.258747122872904`}, {3.`, 5.`}, {4.`,
     0.`}, {4.`, 0.4002986944615309`}, {4.`, 
    0.8326451978805829`}, {4.`, 1.2996052494743657`}, {4.`, 
    1.803950000871885`}, {4.`, 2.348672461377994`}, {4.`, 
    2.9370052598409973`}, {4.`, 3.5724398285307286`}, {4.`, 
    4.258747122872904`}, {4.`, 5.`}, {5.`, 0.`}, {5.`, 
    0.4002986944615309`}, {5.`, 0.8326451978805829`}, {5.`, 
    1.2996052494743657`}, {5.`, 1.803950000871885`}, {5.`, 
    2.348672461377994`}, {5.`, 2.9370052598409973`}, {5.`, 
    3.5724398285307286`}, {5.`, 4.258747122872904`}, {5.`, 5.`}};
melms = {QuadElement[{{1, 11, 12, 2}, {2, 12, 13, 3}, {3, 13, 14, 
      4}, {4, 14, 15, 5}, {5, 15, 16, 6}, {6, 16, 17, 7}, {7, 17, 18, 
      8}, {8, 18, 19, 9}, {9, 19, 20, 10}, {11, 21, 22, 12}, {12, 22, 
      23, 13}, {13, 23, 24, 14}, {14, 24, 25, 15}, {15, 25, 26, 
      16}, {16, 26, 27, 17}, {17, 27, 28, 18}, {18, 28, 29, 19}, {19, 
      29, 30, 20}, {21, 31, 32, 22}, {22, 32, 33, 23}, {23, 33, 34, 
      24}, {24, 34, 35, 25}, {25, 35, 36, 26}, {26, 36, 37, 27}, {27, 
      37, 38, 28}, {28, 38, 39, 29}, {29, 39, 40, 30}, {31, 41, 42, 
      32}, {32, 42, 43, 33}, {33, 43, 44, 34}, {34, 44, 45, 35}, {35, 
      45, 46, 36}, {36, 46, 47, 37}, {37, 47, 48, 38}, {38, 48, 49, 
      39}, {39, 49, 50, 40}, {41, 51, 52, 42}, {42, 52, 53, 43}, {43, 
      53, 54, 44}, {44, 54, 55, 45}, {45, 55, 56, 46}, {46, 56, 57, 
      47}, {47, 57, 58, 48}, {48, 58, 59, 49}, {49, 59, 60, 50}}, {10,
      10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 
     10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 
     10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10}]};
bcEle = {LineElement[{{1, 11}, {2, 1}, {3, 2}, {4, 3}, {5, 4}, {6, 
      5}, {7, 6}, {8, 7}, {9, 8}, {20, 10}, {10, 9}, {11, 21}, {30, 
      20}, {21, 31}, {40, 30}, {31, 41}, {50, 40}, {41, 51}, {51, 
      52}, {52, 53}, {53, 54}, {54, 55}, {55, 56}, {56, 57}, {57, 
      58}, {58, 59}, {59, 60}, {60, 50}}, {2, 1, 1, 1, 1, 1, 1, 1, 1, 
     0, 1, 2, 0, 2, 0, 2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}]};
mesh = ToElementMesh["Coordinates" -> crd, "MeshElements" -> melms, 
   "BoundaryElements" -> bcEle];
mesh["Wireframe"["MeshElement" -> "BoundaryElements", 
  "MeshElementMarkerStyle" -> Blue, 
  "MeshElementStyle" -> {Black, Green, Red}, ImageSize -> Medium]]
mesh["Wireframe"["MeshElementStyle" -> {FaceForm[Red]}, 
  ImageSize -> Medium]]

Mesh

If I apply a modified version of the heat equation taken from Solving PDEs with FEM Tutorial to the mesh, then NDSovleValue says it can't find the Marker:

op = \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[t, x, y]\)\) - \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y}\), \(2\)]\(u[t, x, y]\)\);
dc = DirichletCondition[u[t, x, y] == 1, 
   ElementMarker == bounds["hot"]];
ufunHeat = 
  NDSolveValue[{op == 0, dc, u[0, x, y] == 0}, 
   u, {t, 0, 5}, {x, y} \[Element] mesh];
frames = Table[
   Plot3D[ufunHeat[t, x, y], {x, y} \[Element] mesh, 
    PlotRange -> {0, 1}], {t, 0, 5, 0.5}];
ListAnimate[frames, SaveDefinitions -> True]

If I remove the ElementMarker condition from the DirichletCondition and replace it with a coordinate condition like so

dc = DirichletCondition[u[t, x, y] == 1, y == 0];

Then NDSolve appears to work normally.

I have used this ElementMarker construct on triangular meshes without issue. Should this construct also work with quads? Do I need to specify the marker differently to get NDSolve to recognize it?

$\endgroup$
  • 1
    $\begingroup$ DirchletCondtions use the markers from PointElements not from MeshElements nor from BoundaryElements. If you can not figure it out by tomorrow I'll write an example. IIRC there is one in the ElementMesh generation tutorial. $\endgroup$ – user21 Sep 3 at 19:35
  • 1
    $\begingroup$ In your example you'd need to add PointElement for the dirichlet condition. $\endgroup$ – user21 Sep 3 at 19:37
  • $\begingroup$ Thank you very much. I will try to find that section. That makes sense since the DC is a nodal condition. Perhaps, just union the element incidents of the BoundaryElements? Ultimately, I would like to be able to use the approach to process other meshes. It is much easier for me to detect boundaries by region ID and normals. $\endgroup$ – Tim Laska Sep 3 at 19:45
  • $\begingroup$ Ideally yes, but what should be done a node 2 in this case LineElement[{{1,2},{2,3}},{1,2}]? Any preference? $\endgroup$ – user21 Sep 3 at 19:56
  • $\begingroup$ Do you mean cases where a node is shared by two boundaries? That is a good question. PointElement[Transpose@{{1, 2, 1}}, {1, 2, 3}] allows me to assign multiple markers to node 1, but will ToElementMesh or NDSolve allow it or do I need determine a precedence hierarchy for the nodal BC's? These are good questions that I will need to experiment with. $\endgroup$ – Tim Laska Sep 3 at 20:16
3
$\begingroup$

Here is one way to do it. It is important to realize that DirichletConditions use the markers present in PointElements not the ones in the MeshElements nor the ones in the BoundaryElements. The markers in the BoundaryElements are used exclusively for NeumannVaues and the markers in MeshElements are exlusively used for PDE coefficients that use markers. This is documented in the ElementMesh Generation tutorial.

We start by extracting the PointElements from the mesh you have:

 mesh["PointElements"]
{PointElement[{{1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, \
{10}, {11}, {20}, {21}, {30}, {31}, {40}, {41}, {50}, {51}, {52}, \
{53}, {54}, {55}, {56}, {57}, {58}, {59}, {60}}]}

mesh[
 "Wireframe"["MeshElement" -> "PointElements", 
  "MeshElementIDStyle" -> Black, ImageSize -> Medium]]

enter image description here

Now we add markers to the PointElements

pMarkers = {7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 99, 0, 99, 0, 99, 0, 99, 0, 
99, 0, 0, 0, 0, 0, 0, 0, 0, 0};

pEle = {PointElement[{{1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, \
{10}, {11}, {20}, {21}, {30}, {31}, {40}, {41}, {50}, {51}, {52}, \
{53}, {54}, {55}, {56}, {57}, {58}, {59}, {60}},
pMarkers]};

Recreate the mesh:

mesh = ToElementMesh["Coordinates" -> crd, "MeshElements" -> melms, 
"BoundaryElements" -> bcEle, "PointElements" -> pEle];

We look at the point element markers. Note that the numbering can be completely different from that in the, say, boundary elements.

mesh["Wireframe"["MeshElement" -> "PointElements", 
  "MeshElementMarkerStyle" -> Blue, 
  "MeshElementStyle" -> {Black, Green, Red}, ImageSize -> Medium]]

enter image description here

A different approach is to use the BoundaryMarkerFunction and the PointMarkerFunction documented in the options section of ToBoundaryMesh.

The reason that the point element markers are not populated automatically is the following: Imagine you have a line segment LineElement[{{1,2},{2,3}},{1,2}] which marker should be attributed to the node with index 2? In general markers are a single positive integer for for a point/line/mesh element. So assigning two marker to a single element is currently not an option. If someone has other ideas please let me know. Have a look at the PointMarkerFunction that could be used.

$\endgroup$
  • $\begingroup$ +1 @user21 I didn't know about all information regarding markers in the first paragraph? Is this documented? I think it should be. $\endgroup$ – Pinti Sep 4 at 6:58
  • 1
    $\begingroup$ @Pinti, sure, I added a link to the relevant section. $\endgroup$ – user21 Sep 4 at 7:14
  • $\begingroup$ @user21 Thank you very much for your detailed answer. I makes much sense and was documented, but I did not make the connection. One question I have is there a logic to how the automatic triangle mesh generation creates PointMarkers or is it arbitrary? One FEM tool I use, has a precedence setting that allows a shared PointElement at the intersection of multiple boundaries to take on the highest rank precedence. There is always a good chance that you chose poorly so you can change the precedence without redoing the mesh. $\endgroup$ – Tim Laska Sep 4 at 12:23
  • $\begingroup$ @TimLaska, interesting: how is this rank precedence defined? For the PointMarkers in triangle meshes, I'd need to look that up. I'll get back to you with that. (need to finish something else first) $\endgroup$ – user21 Sep 4 at 12:54
  • $\begingroup$ @TimLaska, so this is done during the call of TriangleLink that provides the actual triangulation. I think this has the same issues we discussed and makes a random choice of point element marker at nodes where it's clear which marker they should belong to. $\endgroup$ – user21 Sep 4 at 14:54

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.