8
$\begingroup$

Previous discussion about this problem has already occurred in this link:

How to get the node id in a specific coordinate in a finite element mesh?

I have the following code that generates a finite element mesh

Needs["NDSolve`FEM`"]
order = 2;
mesh = ToElementMesh[ImplicitRegion[300 <= Sqrt[x x + y y], {x, y}], {{0, 1000}, {0, 1000}}, MaxCellMeasure -> 2000, "MaxBoundaryCellMeasure" -> 10,"MeshOrder" -> order, "NodeReordering" -> True]
mesh["Wireframe"]
topol = mesh["MeshElements"][[1, 1]];
nnodes = mesh["Coordinates"];

But instead to get the the ids of the nodes that are located in the corners, which was the question in the link above, i want all the ids in a specific line in order to impose the boundary conditions to my problem. For example: i need all the node ids of the line located in the curved hole in order to impose an internal presure using Newmamm boudary conditions. How to solve this?

Published work: https://onlinelibrary.wiley.com/doi/full/10.1002/cae.21958

enter image description here

$\endgroup$
  • $\begingroup$ Are you planning to use NDSolve? If so, you don't need node IDs in order to setup boundary conditions. $\endgroup$ – Young Feb 5 '17 at 2:39
  • $\begingroup$ @Young i'm not using NDSolve, i have written my own code. $\endgroup$ – Diogo Feb 5 '17 at 2:46
  • $\begingroup$ reference.wolfram.com/language/FEMDocumentation/tutorial/… $\endgroup$ – Young Feb 5 '17 at 16:06
  • $\begingroup$ @Young i cant use the following functions described in the documentation that you refered : bn = bmesh["BoundaryNormals"]; mean = Mean /@ GetElementCoordinates[bmesh["Coordinates"], #] & /@ ElementIncidents[bmesh["BoundaryElements"]]; Show[ bmesh["Wireframe"], Graphics[ MapThread[ Arrow[{#1, #2}] &, {Join @@ mean, Join @@ (bn/14 + mean)}]]] $\endgroup$ – Diogo Feb 5 '17 at 18:36
  • $\begingroup$ My Mathematica (11) does not recognise "BoundaryNormals" and GetElementCoordinates $\endgroup$ – Diogo Feb 5 '17 at 18:37
8
$\begingroup$

Here is how I'd do it. Boundary conditions are either applied at nodes (Dirichlet type conditions) or over edges/faces (Neumann type conditions). So, I'd generate a mesh that contains markers for points and edges.

For that we first write a helper function that inserts point markers into the mesh:

Needs["NDSolve`FEM`"]
pointMarkerFunction = Compile[{{coords, _Real, 2}},
   Module[{x, y},
      {x, y} = #;
      Which[
       (* order matters *)
       x^2 + y^2 <= 300^2, 5,
       x == 1000., 2,
       x == 0., 4,
       y == 0., 1,
       y == 1000., 3,
       True, 6]] & /@ coords];

And generate a mesh:

mesh = ToElementMesh[
   ImplicitRegion[
    300 <= Sqrt[x x + y y], {x, y}], {{0, 1000}, {0, 1000}}, 
   MaxCellMeasure -> 20000
   , "PointMarkerFunction" -> pointMarkerFunction
   ];

Now we look at the point markers:

Show[mesh["Wireframe"],
 mesh["Wireframe"["MeshElement" -> "PointElements", 
   "MeshElementMarkerStyle" -> Red]]]

enter image description here

There are a few things to note here. First, note that the mid-side nodes of the second order mesh do not yet have the proper markers. That's because the marker value of mid-side nodes is derived from the marker values of the edges, which we have not defined yet. The second thing to note is that the order in the Which statement matters. It's set up in such a way that the markers 5 and 2 include both respective corner points. Depending on which edge the corner points should belong to the order here needs to be adjusted.

Next we add the markers for the edges:

boundaryMarkerFunction = 
  Compile[{{boundaryElementCoords, _Real, 
     3}, {pointMarkres, _Integer, 2}},
   Module[{pm1 = #[[1]], pm2 = #[[2]]},
      Which[
       pm1 == pm2, pm1,
       pm1 == 1 || pm2 == 1, 1,
       pm1 == 3 || pm2 == 3, 3,
       pm1 == 4 || pm2 == 4, 4,
       True, 6 ]] & /@ pointMarkres];

mesh = ToElementMesh[
   ImplicitRegion[
    300 <= Sqrt[x x + y y], {x, y}], {{0, 1000}, {0, 1000}}, 
   MaxCellMeasure -> 20000
   , "BoundaryMarkerFunction" -> boundaryMarkerFunction
   , "PointMarkerFunction" -> pointMarkerFunction
   ];

Show[mesh["Wireframe"],
 mesh["Wireframe"["MeshElement" -> "PointElements", 
   "MeshElementMarkerStyle" -> Red]]]

enter image description here

Now, the mid-side nodes are correct.The boundary marker function gets two arguments; the coordinates of the edge and the point markers of the connected nodes. In this case I choose to use the point markers of the nodes to determine the edge markers.

We can look at the edge markers:

Show[mesh["Wireframe"],
 mesh["Wireframe"["MeshElement" -> "BoundaryElements", 
   "MeshElementMarkerStyle" -> Blue]]]

enter image description here

Here is a function that extracts the coordinates at a specific marker:

getNodes[mesh_, marker_] := Module[{coords, markerpos, inci},
  coords = mesh["Coordinates"];
  markerpos = 
   Flatten[Position[Flatten[ElementMarkers[mesh["PointElements"]]], 
     marker]];
  inci = Flatten[ElementIncidents[mesh["PointElements"]]][[markerpos]];
  coords[[inci]]
  ]

Let's use it for point marker 5:

bcCoords = getNodes[mesh, 5];

And check it:

Norm[300 - Sqrt[Total[bcCoords^2, {2}]]] < 10^-10
True

This gets the curved edges right. As a peak into the future, the marker insertion will be done automatically in a future version of Mathematica. You will find more documentation about markers in the ToBoundaryMesh, ToElementMesh ref. pages and the ElementMesh generation tutorial.

$\endgroup$
4
$\begingroup$

This creates the element mesh from the boundary so that the boundary coordinates are easily extracted. Then using the ID selection from How to get the node id in a specific coordinate in a finite element mesh? in conjunction with the selected boundary coordinates for the arc, I believe gives the desired answer.

Needs["NDSolve`FEM`"]
order = 2;
boundaryMesh = 
 ToBoundaryMesh[
  ImplicitRegion[
   300 <= Sqrt[x x + y y], {x, y}], {{0, 1000}, {0, 1000}},
  "MaxBoundaryCellMeasure" -> 10, "MeshOrder" -> order, 
  "NodeReordering" -> True]
elementMesh = ToElementMesh[boundaryMesh]

elementMesh["Wireframe"]
enodes = elementMesh["Coordinates"];
bnodes = boundaryMesh["Coordinates"];

arcCord = Select[bnodes, #[[1]] <= 300 && #[[2]] <= 300 &];
perCord = Select[bnodes, #[[1]] == 1000 || #[[2]] == 1000 &];
ListPlot[{arcCord, perCord}, AspectRatio -> 1]

findIds[nodes_, coords_] := Nearest[nodes -> Automatic, coords, 1]
arcIds = Flatten@
  Table[findIds[enodes, arcCord[[i]]], {i, 1, Length[arcCord]}]
perIds = Flatten@
  Table[findIds[enodes, perCord[[i]]], {i, 1, Length[perCord]}]
$\endgroup$
  • 2
    $\begingroup$ You already have shown plots which indicate that the points you found are at the correct boundaries. By looking at something like: Show[elementMesh["Wireframe"], Graphics[{Red,Point[enodes[[arcIds]]]}], PlotRange -> {{0, 305}, {0, 305}}] you can also see they are the correct positions. Being able to easily check such things graphically in my opinion is one of the strengths of mathematica which are often underestimated... $\endgroup$ – Albert Retey Feb 5 '17 at 20:32
  • $\begingroup$ I'd be a bit careful: You first generate a boundary mesh and from that a full element mesh. That mesh does not have curved boundaries. For that to happen ToElementMesh will need the ImplicitRegion. In a future version there will be a way to work around this but not in version 11.0 $\endgroup$ – user21 Feb 6 '17 at 8:56
  • $\begingroup$ If you look at Sqrt[Total[enodes[[arcIds]]^2, {2}]] you will see what I mean. $\endgroup$ – user21 Feb 6 '17 at 8:57
  • $\begingroup$ @user21 Is just reducing the boundary cell measure a solution? Can you look at mathematica.stackexchange.com/questions/136372/… and tell me if there is a better way to produce that mesh? I'm primarily concerned with accuracy. $\endgroup$ – Young Feb 6 '17 at 16:19
1
$\begingroup$

Here is my attempt:

First creating the mesh:

Needs["NDSolve`FEM`"]
order = 2;
mesh = ToElementMesh[ImplicitRegion[300 <= Sqrt[x x + y y], {x, y}], {{0,1000}, {0, 1000}}, MaxCellMeasure -> 20000];
nnodes = mesh["Coordinates"];

Then using the method created by @Albert Retey to find the ids, i have found the nodes in the line extremities:

FindIds[nodes_, coords_] :=Block[{nodestofind = Nearest[nodes, coords][[All, 1]]},Flatten[Position[nodes, Alternatives @@ nodestofind]]]
nodetofind = {{300, 0}, {0, 300}};
{hole1, hole2} = FindIds[nnodes, nodetofind]

The following is a method that given two ids in the mesh boundary return all nodes between the specified points:

FindBCNodes[id1_, id2_, nnodes_, order_] := 
  Block[{h1, h2, Flag, A, B, tol, i, norm, linenodes},
   A = nnodes[[id1]];
   B = nnodes[[id2]];
   tol = 10^-3;
   linenodes = {};
   Flag = True;
   h1 = Sqrt[A[[1]] A[[1]] + A[[2]] A[[2]]];
   h2 = Sqrt[B[[1]] B[[1]] + B[[2]] B[[2]]];
   For[i = 1, i <= Length[nnodes], i++,
    (*Vertical Line*)
    If[Abs[A[[1]] - B[[1]]] < tol && 
       Abs[nnodes[[i]][[1]] - A[[1]]] < tol,
      If[Between[
        nnodes[[i]][[2]], {{A[[2]], B[[2]]}, {B[[2]], A[[2]]}}],
       AppendTo[linenodes, nnodes[[i]]];
       ];
      Flag = False;
      ,
      (*Horizontal Line*)
      If[Abs[A[[2]] - B[[2]]] < tol && 
         Abs[nnodes[[i]][[2]] - A[[2]]] < tol,
        If[
         Between[nnodes[[i]][[
           1]], {{A[[1]], B[[1]]}, {B[[1]], A[[1]]}}],
         AppendTo[linenodes, nnodes[[i]]];
         ];
        Flag = False;
        ,
        norm = Norm[nnodes[[i]]];
        If[(Abs[norm - h1]) < tol && Abs[norm - h2] < tol  && 
          Flag != False,
         AppendTo[linenodes, nnodes[[i]]];
          ];
        ];
      ];
    ];
   linenodes
   ];

Calling the function:

bcCoords = FindBCNodes[hole1, hole2, nnodes, order]
Norm[300 - Sqrt[Total[bcCoords^2, {2}]]] < 10^-10
True
$\endgroup$

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.