3D internal subregion
Maybe it helps to have different "RegionMarker"
for different domains? I am not sure if this is useful for OP original problem regarding boundary conditions, but I would like to show another approach to generating such FEM mesh.
Since the geometry of this pipe is quite simple we can mesh it with structured quadrilaterals. This also gives you a lot more control over element size and distribution. I have been developing a package MeshTools to help with this.
Load the package after installing it and define geometry parameters.
Needs["MeshTools`"]
center = {0, 0};
r1 = 0.9; (* inner radius *)
r2 = 1.; (* outer radius *)
n = 10; (* number of elements on quarter of disk perimeter *)
length = 4; (* length of cylinder *)
n2 = 10; (* number of elements in length direction *)
Create structured mesh on Disk
. This represents the "fluid" domain.
diskMesh = DiskMesh[center, r1, n]
(* ElementMesh[{{-0.9, 0.9}, {-0.9, 0.9}}, QuadElement["<" 220 ">"]}] *)
Create structured mesh on Annulus
, which represents the "solid" domain. Element spacing on the edge should match the one of diskMesh
.
annulusMesh = AnnulusMesh[center, {r1, r2}, {4*n, 2}]
(* ElementMesh[{{-1., 1.}, {-1., 1.}}, {QuadElement["<" 80 ">"]}] *)
Observe the result for both meshes (they are still two different ElementMesh
objects).
Show[
diskMesh["Wireframe"["MeshElementStyle" -> FaceForm@Yellow]],
annulusMesh["Wireframe"["MeshElementStyle" -> FaceForm@Orange]]
]
And then the crucial step. Extrude quadrialteral meshes in Z direction, add "RegionMarker" of your choice and merge meshes into one ElementMesh
object.
mesh3D = MergeMesh[{
AddMeshMarkers[ExtrudeMesh[diskMesh, length, n2],"MeshElementsMarker" -> 1],
AddMeshMarkers[ExtrudeMesh[annulusMesh, length, n2],"MeshElementsMarker" -> 2]
}]
mesh3D["Wireframe"[
"MeshElement" -> "MeshElements",
"MeshElementStyle" -> FaceForm /@ {Yellow, Orange},
Axes -> True
]]
I managed to find a dirty solution on my own. If the fluid part is removed completely like what I did in the second try, Mathematica will not consider it as a RegionHole. Thus, the idea is not to remove it completely. Somehow this dirty solution works with the equations too.
wholeRegion[x_, y_, z_] := (x^2 + y^2 <= 0.1^2) && 0 <= z <= 1
fluid[x_,y_,z_] := x^2 + y^2 < (0.9/10)^2 && 0.01 <= z <= 0.99
bmesh = ToBoundaryMesh[
ImplicitRegion[wholeRegion[x, y, z] && ! fluid[x, y, z], {x, y, z}],
"RegionHoles" -> {0, 0, 0.5},
"MaxBoundaryCellMeasure" -> 10^(-4)]
mesh = ToElementMesh[bmesh, "RegionHoles" -> None(*,
MaxCellMeasure\[Rule]10^-5*)]
But still, I hope that Wolfram could optimize the RegionHole to make generating the 3D internal subregion easier.
Update: I have updated the code a bit to work in other cases as well.
Yes, this is too difficult. Here is a workaround for this case. It's not trivial but it works well for the problem at hand. Basically, we generate a boundary mesh from a 2D mesh that has an interior region:
m = ToElementMesh[Annulus[{0, 0}, {9/10, 1}], "RegionHoles" -> None,
"MeshOrder" -> 1, MaxCellMeasure -> Infinity];
m["Wireframe"]
Next, we extract the coordinates and get the largest incident number.
coords = m["Coordinates"];
c1 = Join[coords, ConstantArray[{0.}, {Length[coords]}], 2];
c2 = Join[coords, ConstantArray[{1.}, {Length[coords]}], 2];
ele = m["MeshElements"];
max = Max[ElementIncidents[ele]]
136
We extract the boundary elements from the mesh and construct the sides:
tohigherIndex[type_[inci_, marker___], max_] :=
type[inci + max, marker]
(* old code: not good *)
(* tmpBM = ToBoundaryMesh[m];
sides = QuadElement[Map[Join[#1, Reverse[#1 + max]] &, #]] & /@
ElementIncidents[
Flatten[MeshElementSplitByMarker[tmpBM["BoundaryElements"]]]];*)
sides = QuadElement[Map[Join[#1, Reverse[#1 + max]] &, #]] & /@
ElementIncidents[
Flatten[MeshElementSplitByMarker[m["BoundaryElements"]]]];
We now construct the boundary mesh:
bm = ToBoundaryMesh["Coordinates" -> Join[c1, c2],
"BoundaryElements" ->
Join[ele, tohigherIndex[#, max] & /@ ele, sides]]
bm["Wireframe"]
And generate the full mesh:
ToElementMesh[bm]["Wireframe"]
Hope this helps.
With the updated version I get this for your second mesh: