How to create a surface mesh embedded in 3D?

Question

I am working on a problem where I want to optimize current densities on 2D surfaces embedded in 3D to generate magnetic fields. The current density on a surface can be written as a scalar function multiplied by the curl of the normal vector to the surface.

Ideally, my problem could be solved in the following steps:

1. Create surfaces, probably as boundaries of regions in 3D
2. Generate a mesh on these surfaces
3. Optimize function values on the nodes (I know how to do this step)
4. Create a contour plot of this function on these surfaces (these contours are basically the wire patterns I am aiming for)

My question is about generating the surfaces, meshing them, and plotting the scalar function.

I've tried to come up with a simple example. Consider the following region:

Needs["NDSolve`FEM`"]
rI = 1/2; rA = 1;
reg = RegionProduct[Annulus[{0, 0}, {rI, rA}, {0, Pi/2}], Line[{{0}, {1}}]]

My first problem: The default boundary mesh does not reproduce the region too well, compared to a BoundaryMeshRegion:

bdm = ToBoundaryMesh[reg];
bdr = BoundaryDiscretizeRegion[reg, MaxCellMeasure -> {"Length" -> .15}];
GraphicsRow[{bdm["Wireframe"], bdr}, ImageSize -> Large]

So, I continue with the BoundaryMeshRegion. Since not all surfaces of my initial region reg will carry a current, I manually select one:

points = MeshCells[bdr, 0] /. Point[arg_] :> arg;
lines = MeshCells[bdr, 1] /. Line[arg_] :> arg;
elements = MeshCells[bdr, 2] /. Polygon[arg_] :> arg;
coords = MeshCoordinates[bdr];

radii=coords/.{x_,y_,z_}:>Sqrt[x^2+y^2];
selector[l_List]:=TrueQ[And@@((Abs[#-1.0]<10^-8)&/@radii[[l]])]
selectedPoints=Union[Flatten[Cases[elements,l_List/;selector[l],1]]];
selectedElements=Cases[elements,l_List/;selector[l],1];
selectedCoordinates=coords[[selectedPoints]];
newIndices=Thread[selectedPoints->Range[Length[selectedPoints]]];
(* Reverse is needed to get the incidents in counterclockwise order *)
newElements=Map[Reverse,selectedElements/.newIndices];
newCoords=selectedCoordinates/.{x_Real,y_Real,z_Real}:>{z,ArcTan[x,y]/(Pi/2)};

Plotting everything shows a reasonable structure, with some low quality elements:

Graphics[{
  FaceForm[], EdgeForm[LightRed],
  Polygon[newElements /. i_Integer :> newCoords[[i]]],
  Blue, Text[#, newCoords[[#]]] & /@ Range[Length[newCoords]],
  Black, Text[#, 
     Mean[newElements[[#]] /. i_Integer :> newCoords[[i]]]] & /@ 
   Range[Length[newElements]]
  }, Frame -> True]

Now I convert my manually extracted mesh back to an element mesh:

em1 = ToElementMesh[
  "Coordinates" -> newCoords,
  "MeshElements" -> {TriangleElement[newElements]}
  ];

and generate some function values:

data1 = em1["Coordinates"] /. {z_Real, u_Real} :> Sin[z Pi] Sin[u Pi];
if1 = ElementMeshInterpolation[{em1}, data1, 
  "ExtrapolationHandler" -> {Function[Indeterminate], 
    "WarningMessage" -> False}]

With this poor mesh resolution, the plot is jagged, but it works as expected.

cp1 = ContourPlot[if1[x, y], {x, y} \[Element] em1, Mesh -> None, 
  PlotPoints -> 75, MaxRecursion -> 3, ContourShading -> None, 
  ContourStyle -> Red, Contours -> Range[0.05, 1, 0.1]]

I know how to extract the contours and transform them back to my original surface, but I have no idea how to plot them directly on the surface. SliceContourPlot3D should be able to plot on a region, but I haven't managed to get that to work when function values are defined on the surface only.

Now, my question: Is there an easier way to generate this mesh, ideally with more control over the mesh quality and avoid conversion from BoundaryMeshRegion? How to plot my contours directly on the original surface?

Any insights will be greatly appreciated.

Answer 1

I haven't done the upgrade to 10.4, this is still stuck in our IT department...

Nevertheless, playing with different combinations of ToElementMesh, MeshRegion and RegionBoundary, I've come up with some reasonable results:

For my still rather simple example I use a cylinder with a cut-out.

reg1=Cylinder[{{-2,0,0},{2,0,0}},1];
reg2=Hexahedron[{{-1,-1,0},{-1,1,0},{1,1,0},{1,-1,0},{-1,-1,1},{-1,1,1}, {1,1,1},{1,-1,1}}];
reg=RegionDifference[reg1,reg2];
bdr=BoundaryDiscretizeRegion[reg,Method->Automatic,MaxCellMeasure->{"Length"->0.1}];
Show[bdr,Axes->True,AxesLabel->{"X","Y","Z"},Boxed->True]

As @user21 mentioned, ToBoundaryMeshin V10.4 will be even better, but this is good enough for now. Next step is to generate some dummy values on the mesh:

em=ToElementMesh[bdr,"MeshOrder"->1];
coordinates=em["Coordinates"];
dummy=Quiet[coordinates/.{x_Real,y_Real,z_Real}:>Sqrt[z^2+y^2]Cos[ArcTan[y,z]]Sin[x \[Pi]/2]/.Indeterminate->0.];
if=ElementMeshInterpolation[{em},dummy,"ExtrapolationHandler"->{Function[Indeterminate],"WarningMessage"->False}];

Now, if I plot this InterpolatingFunction over the surface defined with BoundaryDiscretizeRegion, the cut-out surfaces are missing.

But if I transfrom the ElementMesh to a MeshRegion and extract the boundary with RegionBoundary, I get what I expected:

mr=MeshRegion[em];
br=RegionBoundary[mr];
SliceContourPlot3D[if[x,y,z],br,{x,y,z}\[Element]mr,ContourShading->None,Contours->Range[-1.05,1.05,0.1],BoxRatios->{2,1,1}]

Answer 2

Here is a way to get an improved boundary, for this you will need version 12.1 and make use of the new OpenCascadeLink:

Needs["OpenCascadeLink`"]
Needs["NDSolve`FEM`"]
reg1 = Cylinder[{{-2, 0, 0}, {2, 0, 0}}, 1];
reg2 = Cuboid[{-1, -1, 0}, {1, 1, 1}];
reg = RegionDifference[reg1, reg2];
ocr = OpenCascadeShapeBooleanRegion[reg];
bmeshOC = OpenCascadeShapeSurfaceMeshToBoundaryMesh[ocr];
bmeshOC["Wireframe"[
  "MeshElementStyle" -> Directive[FaceForm[Green], EdgeForm[Red]]]]

And then go ahead with:

em = ToElementMesh[bmeshOC, "MeshOrder" -> 1];
coordinates = em["Coordinates"];
dummy = Quiet[
   coordinates /. {x_Real, y_Real, z_Real} :> 
      Sqrt[z^2 + y^2] Cos[ArcTan[y, z]] Sin[x \[Pi]/2] /. 
    Indeterminate -> 0.];
if = ElementMeshInterpolation[{em}, dummy, 
   "ExtrapolationHandler" -> {Function[Indeterminate], 
     "WarningMessage" -> False}];
mr = MeshRegion[em];
SliceContourPlot3D[
 Sqrt[z^2 + y^2] Cos[ArcTan[y, z]] Sin[x \[Pi]/2], mr, 
 Element[{x, y, z}, mr], ContourShading -> None, 
 Contours -> Range[-1.05, 1.05, 0.1], BoxRatios -> {2, 1, 1}, 
 Boxed -> False, Axes -> None]

Answer 3

It seems to me you can get the same result without invoking external packages or using an interpolation function, since you have a surface and you have a function defined at all points in space,

Taking your region bdr,

SliceContourPlot3D[
 Sqrt[z^2 + y^2] Cos[ArcTan[y, z]] Sin[x π/2], bdr, {x, -2, 
  2}, {y, -1, 1}, {z, -1, 1}, ContourShading -> None, 
 Contours -> Range[-1.05, 1.05, 0.1], BoxRatios -> {2, 1, 1}]

gives