I am trying to mesh the following domain to solve a heat transfer + fluid flow problem:

computational domain

The continuity+momentum equations are to be solved in $ABGH$, while the energy equation is to be solved across the entire domain, i.e., $ABCDEFGH$ (basically line $CF$ is the solid-fluid interface). To generate mesh for $ABGH$, I used the following:

L = 0.050, d = 0.003, e = 0.005, delta = 0.010, l = L/d + 2 delta/d
mesh = ToElementMesh[FullRegion[2], {{0, l}, {0, 1}}, MaxCellMeasure -> 8 10^-3];

(*ElementMesh[{{0., 23.3333}, {0., 1.}}, {QuadElement["<" 3132 ">"]}]*)

which produces the following:

ABGH mesh

Now, to mesh the entire domain, I use the tool StructuredMesh in the following code (about which I learned from the various answers on this post):

raster = {{{0, 0}, {delta/d, 0}, {delta/d, -e/d}, {delta/d + L/d, -e/d}, {delta/d + L/d, 0}, {l, 0}}, {{0, 1}, {delta/d, 1}, {delta + L/(3 d), 1}, {delta + (2 L/(3 d)), 1}, {delta/d + L/d, 1}, {l, 1}}}
mesh2 = StructuredMesh[raster, {200, 20}]
(*ElementMesh[{{0., 23.3333}, {-1.66667, 1.}}, {QuadElement["<" 4000 ">"]}]*)

This produces the following:


To produce this, I have basically defined the domain using two lines, each having six points through the raster command. However, as visible the mesh quality is pretty bad.

How should the $ABCDEFGH$ domain be meshed to achieve a smooth quad mesh as $ABGH$? Also, it will be beneficial if I can have the same number of elements in the ABGH portion, while meshing ABCDEFGH as it had when it was individually meshed.

Also, this leads me to think, can both these sub-domains be meshed separetly and then joined?

Update 1 Upon recommendations in the comment, I tried the following:

MergeMesh for different size domains work on BoundaryMesh type object. Atleast, that is what I found in the documentation.

{L = 0.050, d = 0.003, e = 0.005, delta = 0.010, l = L/d + 2 delta/d};
reg1 = ImplicitRegion[0 <= x <= l && 0 <= y <= 1, {x, y}];
reg2 = ImplicitRegion[
   delta/d <= x <= l - delta/d && -e/d <= y <= 0, {x, y}];
reg3 = RegionUnion[reg1, reg2];

bm1 = ToBoundaryMesh[reg1]
bm2 = ToBoundaryMesh[reg2]
bm3 = MergeMesh[{bm1, bm2}]
meshfluid = ToElementMesh[bm1, MaxCellMeasure -> 0.01]["Wireframe"]
meshenergy = ToElementMesh[bm3, MaxCellMeasure -> 0.01]["Wireframe"]
(*meshsolid = ToElementMesh[bm2, MaxCellMeasure -> 0.01]["Wireframe"]*)

This leads to the following:


However, this is a triangular mesh. TriangleToQuadMesh does not work on meshenergy object. Also if you see the image above closely, the refinement at the interface (i.e., internal boundary) is localised to a few regions. Also MeshCellCount does not work when I try it on meshenergy, meshfluid or meshsolid.

Update 2

A comparison among tri and quad mesh results

I used the answer given here by Oleksii, to solve the problem of conjugate heat transfer (reciprocating flow over a heated block). Obviously, I have utlised the domain described in this question. Both the cases have been run for a flow time of 5 seconds, with similar MaxStepSize setting

Triangular mesh results

Created using the Update 1 in this question.

Please note that ToElementMesh creates a 2nd order mesh. To interpolate the p which is supplied with p->1 in NDSolve, I had to create a 1st order mesh using MeshOrderAlteration

  1. Mesh statistics

enter image description here

  1. Solid temperature variation at point ((delta+L)/d/2,-e/d/2)

enter image description here

  1. Fluid temperature variation at both the inlets of the fluid domain

enter image description here

Quad mesh results

Mesh generated using Oleksii's answer in this question

Similar number of elements in y-direction in the fluid domain (compared to trimesh) could be achieved with the setting NyF. Hence, total number of elements is lesser than trimesh.

  1. Mesh statistics

enter image description here

  1. Solid temperature variation at point ((delta+L)/d/2,-e/d/2)

enter image description here

  1. Fluid temperature variation at both the inlets of the fluid domain

enter image description here

In conclusion, the results seem fairly simialr. I will run these tests for longer flow times and report if there are any discrepancies that arise. However, it seems quadmesh allows more control over element size in such regular domains (thus leading to lesser total element count).

Addendum to Update 2 The following contour plots, show the difference in solution between the quad and triangle mesh (at the same flow time). Please note that the triangle mesh is generated using the answer of @user21.

Solid temperature

enter image description here

Fluid temperature Have replaced AspectRatio->Automatic to AspectRatio->1/2, as the fluid region is thin.

enter image description here

Note: Total elements for triangle mesh: 3246 Total elements for quad mesh: 1248

  • 1
    $\begingroup$ You can try to build uniform meshes separately for fluid and solid domains and then to join them by MergeMesh from MeshTools $\endgroup$ Feb 4, 2023 at 12:11
  • $\begingroup$ @OleksiiSemenov I gave it a try. See Update. $\endgroup$
    – Avrana
    Feb 4, 2023 at 13:08
  • 1
    $\begingroup$ Great. A few comments. You should not use the mesh from update 1 - this already looks bad and I suspect it has holes or some such. For the triangular mesh use the mesh that I used for the quad mesh generation. $\endgroup$
    – user21
    Feb 10, 2023 at 8:31
  • 1
    $\begingroup$ If you want to put in the work, what would be much more informative is to plot the difference between the solutions. But again, for this type of problem I would not expect much of a difference in the solution. Triangle elements are typically a bit more efficient to evaluate. $\endgroup$
    – user21
    Feb 10, 2023 at 8:33
  • $\begingroup$ @user21 Thankyou for the comments. See "Addendum to Update 2". The fluid domain shows considerable difference it seems. $\endgroup$
    – Avrana
    Feb 11, 2023 at 3:36

2 Answers 2


Function MergeMesh from MeshTools package allows to join separate meshes. Let's take advantage of this useful function for this problem.

L = 0.050; d = 0.003; e = 0.005; delta = 0.010; 
l = L/d + 2 delta/d;

NxS = 40;
NyS = 10;
MeshSolid = 
  RectangleMesh[{delta/d, -e/d}, {(L + delta)/d, 0}, {NxS, NyS}];

NxF = 10;
NyF = 10;
MeshFluid1 = RectangleMesh[{0, 0}, {delta/d, 1}, {NxF, NyF}];
MeshFluid2 = RectangleMesh[{(delta + L)/d, 0}, {l, 1}, {NxF, NyF}];
MeshFluid3 = RectangleMesh[{(delta)/d, 0}, {(L + delta)/d, 1}, {NxS, NyF}];
MeshFluid = MergeMesh[MergeMesh[MeshFluid1, MeshFluid3], MeshFluid2]

MeshTotal = MergeMesh[MeshFluid, MeshSolid];


  • $\begingroup$ Much much thanks. This is exactly what I wanted. I wasnt aware of the RectangleMesh function too. Keeping similar NxS probably helped. I would also like to ask, if such rectangle meshes are better for fluid flow problems ? $\endgroup$
    – Avrana
    Feb 4, 2023 at 19:20
  • $\begingroup$ @Avrana, it would be interesting to hear what you find when you use this mesh in place of a triangulated region. I doubt it is going to make much of a difference, but I could be wrong. $\endgroup$
    – user21
    Feb 6, 2023 at 14:49
  • 1
    $\begingroup$ @Avrana I will respond on your question about tri mesh vs quads mesh a little bit later. It is blackout in my city now and I am on smartphone $\endgroup$ Feb 6, 2023 at 15:23
  • 1
    $\begingroup$ @Avrana about tri mesh vs quad mesh in cfd. It is popular debatable issue. And as far as I know the choice of proper mesh depends on particular problem in question. May be this work researchgate.net/publication/… will be interesting for you. As concerns the problem you are working on, I think you should carry out the computational experiments with different meshes and then make a choice. Interesting to look at these results. $\endgroup$ Feb 7, 2023 at 21:23
  • 1
    $\begingroup$ @Avrana The points $\{\frac{\delta}{d},0\}$ and $\{\frac{L+\delta}{d},0\}$ are certainly the trouble spots where the large temperature gradients appears. And it would be better to make gradual mesh refinement there. $\endgroup$ Feb 7, 2023 at 21:28

A different approach might be to use ToQuadMesh:

L = 0.050; d = 0.003; e = 0.005; delta = 0.010; l = L/d + 2 delta/d;
bmesh = ToBoundaryMesh[
   "Coordinates" -> {{0, 0}, {delta/d, 
      0}, {delta/d, -e/d}, {delta/d + L/d, -e/d}, {delta/d + L/d, 
      0}, {2*delta/d + L/d, 0}, {L/d + 2*delta/d, 1}, {0, 1}}, 
   "BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4, 
        5}, {5, 6}, {6, 7}, {7, 8}, {8, 1}, {2, 5}}]}];

enter image description here

mesh = ToElementMesh[bmesh];
quadMesh = ToQuadMesh[mesh];

enter image description here

Not, quite what you were looking for but perhaps an alternative option. Note, that the internal boundary is not respected by ToQuadMesh. If you specify region markers then it works:

mesh = ToElementMesh[bmesh, 
   "RegionMarker" -> {{{L/2, 0.9}, 1, 1}, {{L/2, 0.1}, 1, 2}}];
quadMesh = ToQuadMesh[mesh];

enter image description here

  • $\begingroup$ Thankyou for this. I will test this mesh too. I hope the RegionMarker can be used to assign separate properties to each zone? $\endgroup$
    – Avrana
    Feb 10, 2023 at 4:48
  • 2
    $\begingroup$ The RegionMarker can indeed do this. The documentation gives some pretty good examples of use cases. $\endgroup$
    – Dunlop
    Feb 10, 2023 at 5:33

Your Answer

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

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