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..