5
$\begingroup$

I want to create a FEM mesh with an inclusion, but I want to define the coordinates of the edge nodes manually, since I need the nodes for a problem that requires periodic boundary conditions.

enter image description here

The definition of the geometry is not the problem, but I do not know how to incorporate the coordinated of the edge nodes into the mesh.

Is there a possibility to do it with mathematica?

Thanks in advances.

Max

$\endgroup$
3
  • $\begingroup$ Are you trying to find the coordinates of the boundary mesh elements added between the coordinates you specify with the geometry? $\endgroup$
    – tad
    May 11, 2020 at 16:08
  • $\begingroup$ I want to define the nodes on each side of the mesh so that they have on the opposite side a counterpart. Let's say an arbitrary node on the edge has the coordinates (0.0,4.5) then the corresponding node on the opposite edge side should have the coordinates (10.0,4.5), assuming a quadratic mesh geometry with dimensions lets say lx=10 and ly=10. $\endgroup$
    – Max
    May 11, 2020 at 16:40
  • 1
    $\begingroup$ Unlike other FEM codes, Mathematica does not require a one-to-one nodal correspondence for the PeriodicBoundaryCondition. $\endgroup$
    – Tim Laska
    May 11, 2020 at 19:33

2 Answers 2

4
$\begingroup$

I'm not sure if I'm misunderstanding your question, but it sounds like all you want is to control the positions of the nodes at the boundary.

<<NDSolve`FEM`
lpts = Cases[
  Table[{x, y}, {x, -10, 10}, {y, -10, 10}], 
  {a_, b_}/;Abs[a] == 10 || Abs[b] == 10,
  2
];
hexpts = CirclePoints[4, 6];
bmesh = ToBoundaryMesh[
  "Coordinates" -> Join[lpts, hexpts],
  "BoundaryElements" -> {
    LineElement[Partition[Last@FindShortestTour[lpts], 2, 1, {1, 1}]],
    LineElement[Partition[Range[6], 2, 1, {1, 1}] + Length[lpts]]
  },
  "RegionHoles" -> {{0, 0}}
];
mesh = ToElementMesh[bmesh];
Show[
  mesh["Wireframe"],
  Graphics[{
    Red,
    Point@lpts
  }]
]

Mesh with specified edge nodes.

I'm importing the FEM tools from NDSolve as I find they allow more control over mesh creation. I'm generating a matrix of regularly spaced points and then using Cases to select only the boundary points and saving them as lpts. These will be unordered for now. Then I generate the points necessary for the hexagon as hexpts.

ToBoundaryMesh wants a list of all the coordinates under the "Coordinates" argument, so I join lpts and hexpts. Then I add some "BoundaryElements". The first LineElement is the square outer edge. I'm using FindShortestTour to find the ordering that gives me a square and then partitioning this ordering to give me a list of consecutive point numbers like {{1, 2}, {2, 3},...,{n - 1, n}, {n, 1}}. This part wants just the position in the "Coordinates" list and not the actual coordinates themselves. Then I do the same thing for the hexagonal points which will be at the very end of the list, so I add the length of lpts to all of those positions. The last step of creating the boundary mesh is to specify that the hexagon at the centre should be a hole.

I generate a full mesh from this boundary mesh. You should be able to pass this mesh directly to NDSolve as you might any other type of mesh. Finally, I'm displaying the mesh along with the lpts I generated originally to make sure they match the edge nodes. You should be able to customize the lpts to be whatever you like. I tried generating a few different lpts and they always seemed to match the nodes.

$\endgroup$
2
  • $\begingroup$ Thank you very much. Thats basically what I wanted. But there is one more thing... a refinement is is quite easy to incorporate via increment=0.1; lpts = Cases[ Table[{x, y}, {x, -10, 10,increment}, {y, -10, 10, increment}], {a_, b_}/;Abs[a] == 10 || Abs[b] == 10, 2 ]; But a coarsening does not really work out, let's say by increment =2. Do you have any suggestions? $\endgroup$
    – Max
    May 12, 2020 at 10:47
  • $\begingroup$ @Max It's being limited by the default value of MaxCellMeasure. Try changing it to something larger like mesh = ToElementMesh[bmesh, MaxCellMeasure -> 2] or higher and it should work, though there does seem to be an upper limit. If increment is higher than 5, I can't seem to get the mesh to only use those nodes. This is probably due to the hexagon in the centre limiting the maximum cell size. $\endgroup$
    – MassDefect
    May 12, 2020 at 16:19
6
$\begingroup$

Mathematica does not require a one-to-one nodal correspondence for the PeriodicBoundaryCondition to work. However, care must be taken to ensure that PeriodicBoundaryCondition does not share nodes with a DirichletCondition.

Here is an example adapted for an inclusion taken from the documentation for PeriodicBoundaryCondition. Please note that there can be artefacts introduced due to implied NeumannConditions on the "source" boundary as discussed in this MSE post. That is why I applied forward and reverse PBC's. It seemed to work.

Needs["NDSolve`FEM`"]
{length, height, xc, yc, r} = {1, 2, 0, 0, 1/8};
{sx, sy, fx, fy} = {-length/2, -height/2, length/2, height/2};
disk = Region@Disk[{xc, yc}, r];
Ω = 
  RegionDifference[Rectangle[{sx, sy}, {fx, fy}], disk];
mesh = ToElementMesh[Ω, MaxCellMeasure -> 0.0005, 
   AccuracyGoal -> 5];
pde = ((Inactive[
       Div][(-{{1, 0}, {0, 1}}.Inactive[Grad][u[x, y], {x, y}]), {x, 
       y}]) - If[1/4 fx <= x <= 3/4 fx && sy/4 <= y <= fy/4, 1., 
     0.] == 0)
Subscript[Γ, D] = 
  DirichletCondition[
   u[x, y] == 0, (y <= sy || y >= fy) && sx < x <= fx];
pbcf = PeriodicBoundaryCondition[u[x, y], x == sx && sy <= y <= fy, 
   TranslationTransform[{length, 0}]];
pbcr = PeriodicBoundaryCondition[u[x, y], x == fx && sy <= y <= fy, 
   TranslationTransform[{-length, 0}]];
ufun = NDSolveValue[{pde, pbcf, pbcr, Subscript[Γ, D]}, 
   u, {x, y} ∈ mesh];
cp = ContourPlot[ufun[x, y], {x, y} ∈ mesh, 
  ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]
Show[MapAt[Translate[#, {length, 0}] &, cp, 1], cp, 
 MapAt[Translate[#, {-length, 0}] &, cp, 1], PlotRange -> All]

Results

For completeness, I show the artefact with specifying only one PBC resulting in a no flux condition on the source wall.

pbc = PeriodicBoundaryCondition[u[x, y], x == sx && sy <= y <= fy, 
   TranslationTransform[{length, 0}]];
ufun = NDSolveValue[{pde, pbc, Subscript[Γ, D]}, 
   u, {x, y} ∈ mesh];
cp = ContourPlot[ufun[x, y], {x, y} ∈ mesh, 
  ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]
Show[MapAt[Translate[#, {length, 0}] &, cp, 1], cp, 
 MapAt[Translate[#, {-length, 0}] &, cp, 1], PlotRange -> All]

Artefact

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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