3
$\begingroup$

I write a code that generates the following mesh. I would like to find out the elements that completely lie within a certain circular region inside this meshed area. How do I do that?

enter image description here

The code that I wrote is shown below:

w=100; (* Box width *)
h=100; (* Box height *)
xmin = 0; 
ymin = 0; 
xmax = xmin + w; 
ymax = ymin + h;
d = 10;(* Mesh Size *)

<< NDSolve`FEM`
(* Create the coordinates *)
coords = {};(* Nodal Coordinate Matrix *)
nx = w/d + 1;(* Number of nodes in the x-direction *)
ny = h/d + 1;(* Number of nodes in the y-direction *)
For[i = 0, 
 i <= nx - 1, i++,
 For[j = 0, j <= ny - 1, j++,
  coords = Append[coords, {xmin + i*d, ymin + j*d}]; 
  ]
 ]

(* Create Element Connectivity Matrix to generate triangular elements *)
elMat = {};(* Element Connectivity Matrix *)
For[i = 1, 
 i <= ny - 1, i++,
 For[j = 1, j <= nx - 1, j++,
  elMat = 
   Append[elMat, {(i - 1)*ny + j, (i - 1)*ny + j + 1, i*ny + j}];
  elMat = Append[elMat, {i*ny + j, i*ny + j + 1, (i - 1)*ny + j + 1}]
  ]
 ]

(* Create Mesh *)
mesh = ToElementMesh["Coordinates" -> coords, 
   "MeshElements" -> {TriangleElement[elMat]}];
mesh["Wireframe"]
$\endgroup$
6
$\begingroup$

The following looks for coordinates the lie within 20 unit radius of the point {50,50}:

elmcrd = GetElementCoordinates[mesh["Coordinates"], #] & /@ 
   ElementIncidents[mesh["MeshElements"]];
mf = RegionMember[Disk[{50, 50}, 20]];
caf = ContainsAll[{True, True, True}];
PositionIndex[caf /@ Flatten[mf /@ elmcrd, 1]]
(* <|False -> {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 
   15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
    32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 
   48, 49, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,
    75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 97, 98, 99, 100, 101, 
   102, 103, 104, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 
   126, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 
   146, 147, 148, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 
   162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 
   175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 
   188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200}, 
 True -> {50, 51, 52, 68, 69, 70, 71, 72, 73, 74, 86, 87, 88, 89, 90, 
   91, 92, 93, 94, 95, 96, 105, 106, 107, 108, 109, 110, 111, 112, 
   113, 114, 115, 127, 128, 129, 130, 131, 132, 133, 149, 150, 151}|> *)

You should note that the ordering of your triangles gives you a negative mesh quality. You should consider re-ordering the triangle elements.

$\endgroup$
5
$\begingroup$

You can do it easily with MeshTools package.

Get["MeshTools`"]

mesh = QuadToTriangleMesh@RectangleMesh[{0, 0}, {100, 100}, {10, 10}]
mesh["Wireframe"]

mesh_rectangle

reg = Disk[{50, 50}, 40];
selected = SelectElements[mesh, RegionMember[reg][{#1, #2}] &]

Show[
 selected["Wireframe"],
 Graphics[{EdgeForm[Red], FaceForm[], reg}]
]

mesh_disk

$\endgroup$
3
  • 1
    $\begingroup$ good to see you ;-) $\endgroup$ – user21 May 8 '20 at 7:25
  • $\begingroup$ This works very well. I have an additional question: Right now, each quad element is split into 2 triangles. Is it possible to split each quad element into 4 triangular elements? Put another way, is it possible to split each quad element using both its diagonals? $\endgroup$ – Pinkesh May 10 '20 at 19:02
  • $\begingroup$ "MeshTools" pacakge doesn't have such option implemented. But of course, splitting triangles again is possible, but you have to write the function yourself. $\endgroup$ – Pinti May 11 '20 at 6:26

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.