# Mesh unit square with triangles?

**URL:** <https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271>\
**Category:** Questions\
**Tags:** graphics-and-visualization, numerical-computation, packages-and-paclets\
**Created:** [September 21, 2016, 3:47am UTC](https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271 "2016-09-21T03:47:35Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![scinar](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/scinar/32/39519_2.png) [@scinar](https://community.wolfram.com/u/scinar)\
**Post date:** [September 21, 2016, 3:47am UTC](https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271/1 "2016-09-21T03:47:35Z")

</div>

Hi All, I am trying to mesh unit square with triangles as in the picture (Later I’ll divide pts by n to get unit square ). TriangleElement can be made by hand if n is small say 2 or 4. But in my case n=64. How can I make a list for TriangleElement? Any suggestion, thanks in advance.

```
n = 1;
    pts = Flatten[Table[{i, j}, {j, 0, n}, {i, 0, n}], 1]
ListPlot[Labeled[#, #] & /@ pts, Axes -> False]
Needs["NDSolve`FEM`"]
    mesh = ToElementMesh["Coordinates" -> pts, "MeshElements" -> {TriangleElement[{{1, 2, 4}, {4, 3, 1}}]}]
    mesh["Wireframe"]

```

---

<div class="post-metadata">

**Author:** ![henrikschachner](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/henrikschachner/32/32119_2.png) [@henrikschachner](https://community.wolfram.com/u/henrikschachner)\
**Post date:** [September 21, 2016, 12:25pm UTC](https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271/2 "2016-09-21T12:25:29Z")

</div>

Hi, here is a very straight forward approach:

```
n = 10; (* to be choosen *)
pts = Flatten[Table[{x, y}, {y, 0, n}, {x, 0, n}], 1];
indexPts = Flatten[Table[{x, y}, {y, 0, n + 1}, {x, 0, n + 1}], 1];

mkTriangs[pt : {x_, y_}] := 
 Sequence[{pt, pt + {0, 1}, pt + {1, 1}}, {pt, pt + {1, 0}, pt + {1, 1}}]

triangs = mkTriangs /@ pts;
indxRules = MapIndexed[Rule[#1, First[#2]] &, indexPts];
index = triangs /. indxRules;

Needs["NDSolve`FEM`"]
elementMesh = 
  ToElementMesh["Coordinates" -> indexPts, "MeshElements" -> {TriangleElement[index]}] // Quiet;

elementMesh["Wireframe"]

```

which results in:

 ![elementMesh.png](https://global.discourse-cdn.com/wolfram/original/3X/4/e/4ece14498e6d1adb230992b0eb4978b9be33f1cd.png)

Regards – Henrik  
[1]: [http://community.wolfram.com//c/portal/getImageAttachment?filename=elementMesh.png&userId=32203](http://community.wolfram.com//c/portal/getImageAttachment?filename=elementMesh.png&userId=32203)

---

<div class="post-metadata">

**Author:** ![sghuisman](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/sghuisman/32/27850_2.png) [@sghuisman](https://community.wolfram.com/u/sghuisman)\
**Post date:** [September 21, 2016, 1:33pm UTC](https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271/3 "2016-09-21T13:33:35Z")

</div>

Needs[“NDSolve`FEM`”]  
ClearAll[MakeMesh]  
MakeMesh[n\_Integer?Positive]:=Module[{pts,tr1,tr2},  
pts=Join@@Table[{i,j},{j,0,n},{i,0,n}];  
tr1=Join@@Table[1+(n+1) j+i+{0,1,2+n},{j,0,n-1},{i,0,n-1}];  
tr2=Join@@Table[1+(n+1) j+i+{0,2+n,1+n},{j,0,n-1},{i,0,n-1}];  
ToElementMesh[“Coordinates”-\>pts,“MeshElements”-\>{TriangleElement[tr1~Join~tr2]}]  
]  
MakeMesh[6][“Wireframe”]

 ![ScreenShot2016-09-21at15.33.22.png](https://global.discourse-cdn.com/wolfram/original/3X/9/d/9dc6c21f84b9a567be7d2b78163a9d2f349e4233.png)

---

<div class="post-metadata">

**Author:** ![sghuisman](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/sghuisman/32/27850_2.png) [@sghuisman](https://community.wolfram.com/u/sghuisman)\
**Post date:** [September 21, 2016, 1:34pm UTC](https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271/4 "2016-09-21T13:34:14Z")

</div>

i didn’t see your reply, before I posted mine, sorry about that. Very similar code!

---

<div class="post-metadata">

**Author:** ![scinar](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/scinar/32/39519_2.png) [@scinar](https://community.wolfram.com/u/scinar)\
**Post date:** [September 21, 2016, 4:14pm UTC](https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271/5 "2016-09-21T16:14:50Z")

</div>

Wow, Awesome!! Thank you so much both of you

---

<div class="post-metadata">

**Author:** ![scinar](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/scinar/32/39519_2.png) [@scinar](https://community.wolfram.com/u/scinar)\
**Post date:** [September 22, 2016, 2:23am UTC](https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271/6 "2016-09-22T02:23:30Z")

</div>

We have a problem here. My ultimate goal is to solve Poisson equation

-\Delta u=f,\quad(x,y)\in\Omega=(0,1)\times(0,1)

u=0,\quad \Gamma=\partial\Omega

where f(x,y)=2\pi^2\sin(\pi x)\sin(\pi y) using FEM and grid as above. (Exact solution is u(x,y)=\sin(\pi x)\sin(\pi y ))

Let’s look at Mathematica’s example..

```
op = -Laplacian[u[x, y], {x, y}] - 20
\[CapitalGamma] = DirichletCondition[u[x, y] == 0, x == 0 || x == 2]
 mesh = ToElementMesh["Coordinates" -> {{0., 0.}, {1., 0.}, {2., 0.}, {2., 1.}, {1.,1.}, {0., 1.}}, "MeshElements" -> {TriangleElement[{{1, 2, 5}, {5, 6, 1}, {2, 3,4}, {4, 5, 2}}]}]
ufun = NDSolveValue[{op == 0, \[CapitalGamma]}, u[x, y], {x, y} \[Element] mesh];
Show[ContourPlot[ufun, {x, y} \[Element] mesh, AspectRatio -> Automatic],
 mesh["Wireframe"]]

```

it solves the problem..

But when I make the same mesh using below code, it is not working

```
n = 2;(*to be choosen*)
pts = Flatten[Table[{x, y}, {y, 0, n - 2}, {x, 0, n - 1}], 1];
indexPts = Flatten[Table[{x, y}, {y, 0, n - 1}, {x, 0, n}], 1];
mkTriangs[pt : {x_, y_}] := 
 Sequence[{pt, pt + {0, 1}, pt + {1, 1}}, {pt, pt + {1, 0}, 
   pt + {1, 1}}]
triangs = mkTriangs /@ pts;
indxRules = MapIndexed[Rule[#1, First[#2]] &, indexPts];
index = triangs /. indxRules;
index = triangs /. indxRules;
Needs["NDSolve`FEM`"]
mesh = ToElementMesh["Coordinates" -> indexPts, 
    "MeshElements" -> {TriangleElement[index]}] // Quiet;
mesh["Wireframe"]

```

Mathematica uses outside first, I don’t know it is matter. If it is matter how can I get this order..

```
ToElementMesh["Coordinates" -> {{0., 0.}, {1., 0.}, {2., 0.}, {2., 1.}, {1.,1.}, {0., 1.}}, "MeshElements" -> {TriangleElement[{{1, 2, 5}, {5, 6, 1}, {2, 3,4}, {4, 5, 2}}]}]

```

Any suggestion..

---

<div class="post-metadata">

**Author:** ![henrikschachner](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/henrikschachner/32/32119_2.png) [@henrikschachner](https://community.wolfram.com/u/henrikschachner)\
**Post date:** [September 23, 2016, 12:00pm UTC](https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271/7 "2016-09-23T12:00:54Z")

</div>

Hi Okkes,

according to my (poor!) understanding here, the “correct ordering” of triangle points means that they should make a left turn. So I simply changed my function `mkTriangs[]`, which now reads:

```
mkTriangs[pt : {x_, y_}] := 
 Sequence[{pt, pt + {1, 0}, pt + {1, 1}}, {pt, pt + {1, 1}, pt + {0, 1}}]

```

with this it seems to work.

Regards – Henrik

---

<div class="post-metadata">

**Author:** ![scinar](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/scinar/32/39519_2.png) [@scinar](https://community.wolfram.com/u/scinar)\
**Post date:** [September 24, 2016, 3:15am UTC](https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271/8 "2016-09-24T03:15:46Z")

</div>

Yes, order of vertex coordinate of each triangle must be in counterclockwise. Thanks for your helps. Here is my final result. Now how can I extract nSol result for some values from NDSolveValue.? I mean instead of graphing it I need data to find error. By the way is there a good way to use Module function to calculate nSol for different n values. n=8,16,32,64. Thanks.

```
Needs["NDSolve`FEM`"]

n = 8;
pts = Join @@ Table[{i, j}, {j, 0, 1, 1/n}, {i, 0, 1, 1/n}];
upperLeft =1/n Flatten[Drop[Table[{{i, j}, {i + 1, j + 1}, {i, j + 1}}, {j, 0, n}, {i, 0,n}], -1, -1], 1];
lowerRight = 1/n Flatten[Drop[Table[{{i, j}, {i + 1, j}, {i + 1, j + 1}}, {j, 0, n}, {i, 0,n}], -1, -1], 1];
triangs = Riffle[upperLeft, lowerRight];
indxRules = MapIndexed[Rule[#1, First[#2]] &, pts];
element = triangs /. indxRules;
mesh = ToElementMesh["Coordinates" -> pts, "MeshElements" -> {TriangleElement[element]}] // Quiet;
mesh["Wireframe"];

op = -Laplacian[u[x, y], {x, y}] - 2 \[Pi]^2 Sin[\[Pi] x] Sin[\[Pi] y];
\[CapitalGamma] = DirichletCondition[u[x, y] == 0., True];
nSol = NDSolveValue[{op == 0., \[CapitalGamma]}, u[x, y], {x, y} \[Element] mesh];

Plot3D[nSol, {x, 0, 1}, {y, 0, 1}]

ContourPlot[nSol, {x, y} \[Element] mesh, AspectRatio -> Automatic]

Show[ContourPlot[nSol, {x, y} \[Element] mesh, AspectRatio -> Automatic], mesh["Wireframe"]]

```

---

<div class="post-metadata">

**Author:** ![scinar](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/scinar/32/39519_2.png) [@scinar](https://community.wolfram.com/u/scinar)\
**Post date:** [September 24, 2016, 3:22am UTC](https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271/9 "2016-09-24T03:22:54Z")

</div>

I would like to change title to Solving Possion Equation using FEM on unit square with right triangle mesh

---

<div class="post-metadata">

**Author:** ![scinar](https://yyz1.discourse-cdn.com/wolfram/user_avatar/community.wolfram.com/scinar/32/39519_2.png) [@scinar](https://community.wolfram.com/u/scinar)\
**Post date:** [September 24, 2016, 10:43pm UTC](https://community.wolfram.com/t/mesh-unit-square-with-triangles/9271/10 "2016-09-24T22:43:17Z")

</div>

Here how we can extract data from NDSolveValue

```
exact[x_, y_] := Sin[\[Pi] x] Sin[\[Pi] y]
exactData = Partition[{{#1, #2}, exact[#1, #2]} & @@@ pts // Flatten, 3] // N
numData = Partition[NDSolveValue[{op ==0., \[CapitalGamma]}, {{#1, #2}, u[#1, #2]} & @@@ pts, {x, y} \[Element] mesh]//Flaten, 3] // N
ListPointPlot3D[{exactData, numData }]

```