Solución:
Basado en el tutorial aquí
Las mallas QuadElement se comportan exactamente igual que las mallas TriangleElement, con la excepción de que, para elementos quad lineales, se necesitan cuatro incidentes por elemento y, para elementos cuadráticos, se necesitan ocho incidentes por elemento.
aquí hay otra posible respuesta:
Needs["NDSolve`FEM`"]
rectangleA = {{0, 4}, {8, 8}};
rectangleB = {{5, 0}, {8, 4}};
pitch = 0.25;
rectAx = 1/pitch (rectangleA[[2, 1]] - rectangleA[[1, 1]]) + 1;
rectAy = 1/pitch (rectangleA[[2, 2]] - rectangleA[[1, 2]]) + 1;
rectBx = 1/pitch (rectangleB[[2, 1]] - rectangleB[[1, 1]]) + 1;
rectBy = 1/pitch (rectangleB[[2, 2]] - rectangleB[[1, 2]]) + 1;
coordinatesA =
Flatten[Table[{x, y}, {x, rectangleA[[1, 1]], rectangleA[[2, 1]], pitch},
{y, rectangleA[[1, 2]], rectangleA[[2, 2]], pitch}], 1];
coordinatesB =
Flatten[Table[{x, y}, {x, rectangleB[[1, 1]], rectangleB[[2, 1]], pitch},
{y, rectangleB[[1, 2]], rectangleB[[2, 2]], pitch}], 1];
coordinates = Join[coordinatesA, coordinatesB];
incidents = Join[
Flatten[
Table[{j*rectAy + i, j*rectAy + i + 1, (j - 1)*rectAy + i + 1, (j - 1)*rectAy + i},
{i, 1, rectAy - 1},
{j, 1, rectAx - 1}], 1],
Flatten[
Table[{j*rectBy + i, j*rectBy + i + 1, (j - 1)*rectBy + i + 1, (j - 1)*rectBy + i},
{i, 1 + Length[coordinatesA], Length[coordinatesA] + rectBy - 1},
{j, 1, rectBx - 1}], 1]
];
mesh1 = ToElementMesh["Coordinates" -> coordinates,
"MeshElements" -> {QuadElement[incidents]},
"NodeReordering" -> True];
mesh2 = MeshOrderAlteration[mesh1, 2];
mesh2["MeshOrder"]
Show[
DiscretizeGraphics[
GraphicsComplex[{{0, 4}, {5, 4}, {5, 0}, {8, 0}, {8, 8}, {0, 8}},
Polygon[{1, 2, 3, 4, 5, 6}]]],
mesh2["Wireframe"],
ListPlot[coordinates]
]
Problema de muestra:
sol = NDSolveValue[{D[u[x, y], x, x] + D[u[x, y], y, y] == 0,
DirichletCondition[u[x, y] == 100, y == 4 && x <= 5 || y <= 4 && x == 5],
DirichletCondition[u[x, y] == 0, y == 8 || x == 8]},
u, {x, y} ∈ mesh2]
sol["ElementMesh"]
ContourPlot[sol[x, y], {x, y} ∈ mesh2, Mesh -> None,
ColorFunction -> "Rainbow", PlotRange -> All,
PlotLegends -> Automatic]
Aquí hay un enfoque diferente. La idea es generar primero una malla triangular de segundo orden. A continuación, se agrega una coordenada central en cada triángulo. Luego dividimos cada triángulo en tres quads de primer orden haciendo uso de la coordenada central recién agregada y los incidentes de triángulo inicial de segundo orden.
Entonces agregamos un nodo en el centro del elemento y creamos los elementos cuádruples {c, 6,1,4}, {c, 4,2,5} y {c, 5,3,6}.
Comience con una malla triangular inicial:
Needs["NDSolve`FEM`"]
m = ToElementMesh[
DiscretizeGraphics[
GraphicsComplex[{{0, 4}, {5, 4}, {5, 0}, {8, 0}, {8, 8}, {0, 8}},
Polygon[{1, 2, 3, 4, 5, 6}]]]];
Dos funciones auxiliares que funcionan para elementos de malla triangular de segundo orden. Si uno desea extender esto a 3D, deberá agregar una función tetToHex.
getMidEleCoords = Compile[{{eleCoords, _Real, 2}},
Mean[eleCoords], RuntimeAttributes -> {Listable}];
triToQuad = Compile[{{triInci, _Integer, 1}, {midInci, _Integer, 0}}, {
{midInci, triInci[[6]], triInci[[1]], triInci[[4]]},
{midInci, triInci[[4]], triInci[[2]], triInci[[5]]},
{midInci, triInci[[5]], triInci[[3]], triInci[[6]]}
}
, RuntimeAttributes -> {Listable}
];
getConverterCode[TriangleElement] := triToQuad;
getNumOfResEle[TriangleElement] := 3
getConvertedEleType[TriangleElement] := QuadElement
Primero, calculamos las nuevas coordenadas dentro de los elementos triángulo / tet:
coords = m["Coordinates"];
numOfCoords = Length[coords];
meshEle = m["MeshElements"];
inci = ElementIncidents[meshEle];
markers = ElementMarkers[meshEle];
eleCoords = GetElementCoordinates[coords, #] & /@ inci;
midEleCoords = getMidEleCoords[eleCoords];
numOfMidCoords = Length /@ midEleCoords;
temp = Join[{numOfCoords}, numOfMidCoords];
temp = FoldList[Plus, temp] + 1;
midEleInci = MapThread[Range, {Most[temp], Rest[temp - 1]}];
newCoords = Join[coords, Join @@ midEleCoords];
(*Max[midEleInci]===Length[newCoords]*)
Luego generamos los nuevos elementos:
newEle = Table[
Block[{thisMeshEle, thisType, thisInci, thisMarker, converter,
newEle, newEleMarkers},
thisMeshEle = meshEle[[i]];
thisType = Head[thisMeshEle];
thisInci = ElementIncidents[thisMeshEle];
thisMarker = ElementMarkers[thisMeshEle];
converter = getConverterCode[thisType];
newEle = Join @@ converter[thisInci, midEleInci[[i]]];
newEleMarkers =
Join @@ Transpose[
ConstantArray[thisMarker, getNumOfResEle[thisType]]];
getConvertedEleType[thisType][newEle, newEleMarkers]
]
, {i, Length[meshEle]}];
Genere una malla cuádruple:
mq1 = ToElementMesh["Coordinates" -> newCoords,
"MeshElements" -> newEle]
(* ElementMesh[{{0., 8.}, {0., 8.}}, {QuadElement["<" 1338 ">"]}] *)
Mira la malla:
mq1["Wireframe"]
Esta es ahora una malla de primer orden:
mq1["MeshOrder"]
1
Uno puede usar MeshOrderAlteration[mq1, 2] para luego generar posteriormente una malla de segundo orden.
Histogram[Join @@ mq1["Quality"]]
La calidad se ve bien, pero valdría la pena probar algún tipo de suavizado (por ejemplo, suavizado laplaciano).
(Actualizar: Hay una pregunta de seguimiento que analiza la idea del suavizado laplaciano. Spoiler: no parece que valga la pena)
Resuelva un problema de modelo sobre la malla:
sol = NDSolveValue[{-Laplacian[u[x, y], {x, y}] == 1,
DirichletCondition[u[x, y] == 0, True]}, u, {x, y} [Element] mq1];
Plot3D[sol[x, y], {x, y} [Element] sol["ElementMesh"]]
Dado que esto generará tres elementos cuádruples por cada triángulo, la medida de la celda de la malla es menor de lo que debe ser, se podría ajustar haciendo uso de “MaxCellMeasure” en la malla del triángulo inicial.
Espero que esto ayude.
En general, eso no es posible en V11.0. La única primitiva a partir de la cual se puede generar una malla cuádruple es
ToElementMesh[Rectangle[]]
o
ToElementMesh[FullRegion[2], {{0, 1}, {0, 1}}]
Entonces, para obtener una malla cuádruple para una forma más complicada, necesitaría programar un poco y usar
ToElementMesh["Coordinates"->..., "MeshElements"->{QuadElement[....]}]
Como otro enfoque, podría fusionar 2 triángulos en un cuádruple y así obtener una malla en su mayoría cuádruple. Si haces eso. Me interesaría ver el resultado.
Aquí hay otra información: si tiene una malla de primer orden, puede hacer esa segunda orden.
m = ToElementMesh["Coordinates" -> coordinates,
"MeshElements" -> {QuadElement[incidents]}];
m["MeshOrder"]
1
m2 = MeshOrderAlteration[m, 2];
m2["MeshOrder"]
2
Dado que su caso no implica límites curvos, no es necesario mover los nodos del lado medio (nodos de segundo orden) a la posición adecuada.