Esta inquietud se puede abordar de diferentes formas, sin embargo te mostramos la solución más completa en nuestra opinión.
Solución:
Demasiado tiempo para un comentario:
Bien, primero configure su sistema para que solo tenga BC no periódicos. Luego mire el tutorial de programación de elementos finitos y use NDSolve ProcessEquations y siga los pasos hasta la llamada a DiscretizePDE y DiscretizeBoundaryConditions. En este punto puede extraer las matrices del sistema. Implemente las condiciones de contorno (no periódicas). Hasta aquí todo está documentado.
Las condiciones de contorno periódicas vienen en pares; una restricción algebraica que agregó al sistema de ecuaciones.
Aquí hay un ejemplo: Cree una matriz de sistema y un vector de carga:
n = 5;
matrix = Table[
Switch[ j - i, -1, 1, 0, -2, 1, 1, _, 0], i, n, j, n ];
loadVector = Table[ i, i, n ];
Ahora, por ejemplo, nos gustaría cambiar las ecuaciones del sistema de modo que u1 == u3. Hacemos esto con multiplicadores lagrangianos:
lv = ConstantArray[0, n];
(* set u1-u3, for example *)
lv[[1]] = 1;
lv[[3]] = -1;
(* add lagrange row *)
matrix2 = Join[matrix, lv];
(* add lagrange col *)
matrix3 = Join[matrix2, Partition[Join[lv, 0], 1], 2]
(*
-2, 1, 0, 0, 0, 1, 1, -2, 1, 0, 0, 0, 0, 1, -2, 1, 0, -1, 0,
0, 1, -2, 1, 0, 0, 0, 0, 1, -2, 0, 1, 0, -1, 0, 0, 0
*)
El sistema ahora es más grande porque agregamos la ecuación u1 – u3 y también necesitamos ampliar el vector de carga; queremos la ecuación u1 – u3 == 0
(* set u1-u3 == 0 *)
load2 = Join[loadVector, 0];
Resuelve las ecuaciones:
res = LinearSolve[matrix3, load2]
-(31/4), -(35/4), -(31/4), -(19/2), -(29/4), -(23/4)
Y comprueba que u1 == u3
res[[1, 3]]
-(31/4), -(31/4)
Espero que esto ayude.
Bien, después de la actualización, así es como lo haría:
reg = Rectangle[0, 0, 1, 1];
(*Obtain FEM state*)
state =
NDSolve`ProcessEquations[eq0 == 0, conds, u,
Element[x1, x2, reg],
Method -> "FiniteElement",
"MeshOptions" -> "MaxCellMeasure" -> 0.2, "MeshOrder" -> 1,
"IntegrationOrder" -> 2];
femstate = state["FiniteElementData"];
methodData = femstate["FEMMethodData"];
initBCs = femstate["BoundaryConditionData"];
initCoeff = femstate["PDECoefficientData"];
sd = state["SolutionData"][[1]];
discretePDE = DiscretizePDE[initCoeff, methodData, sd];
stiffness = discretePDE["StiffnessMatrix"];
load = discretePDE["LoadVector"];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];
DeployBoundaryConditions[load, stiffness, discreteBCs]
Hasta aquí todo bien. El ajuste “InterpolationOrder” que utilizó es para poder especificar un orden de interpolación diferente para sistemas acoplados, que no tenemos aquí. Entonces, el orden de interpolación siempre será el mismo que el orden de la malla. Además, extraemos todo de los datos estatales. Eso puede ser un poco más conveniente que configurarlo manualmente. Pero, por supuesto, haz lo que quieras.
Extrae y visualiza la malla:
mesh = NDSolve`SolutionDataComponent[sd, "Space"]["ElementMesh"];
Show[
mesh["Wireframe"["MeshElementIDStyle" -> Blue]]
, mesh["Wireframe"["MeshElement" -> "PointElements",
"MeshElementIDStyle" -> Red]]]
Ahora, extraemos las posiciones de las coordenadas relevantes:
c = mesh["Coordinates"];
cornerPos = Position[c, Alternatives @@ Tuples[0., 1., 2]];
tc = Transpose[c];
zeroPos = Position[tc[[1]], 0.]
onePos = Position[tc[[1]], 1.]
(*
1, 2, 3, 4
13, 14, 15, 16
*)
left = Flatten[Complement[zeroPos, cornerPos]]
right = Flatten[Complement[onePos, cornerPos]]
(*
2, 3
14, 15
*)
Usamos el left nodos como nodos maestros y el right nodos como nodos esclavos.
dof = methodData["DegreesOfFreedom"];
lmdof = Length[left];
Ahora, vamos a hackear un poco. Esto no está documentado. Vamos a crear un DiscretizedBoundaryConditionData que luego podemos usar en DeployBoundaryConditions. Esto es un poco arriesgado ya que la estructura de datos internos puede cambiar en el futuro. (Es posible que desee envolver esto en una función y si la estructura de datos cambia, entonces solo necesitará arreglarlo en esa función y no en todos los lugares donde usó el código).
Creamos una matriz donde las posiciones en cada fila se establecen en +1 y -1 como se explicó anteriormente.
diriMat = SparseArray[, lmdof, dof];
MapThread[(diriMat[[#1, #2]] = 1, -1) &, Range[lmdof],
Transpose[left, right]];
Especifique dónde tenemos las filas de Dirichlet:
diriRows = left;
Especifique los valores de Dirichlet que son cero en nuestro caso:
diriVals = SparseArray[, lmdof, 1];
Creamos el DiscretizedBoundaryConditionData estructura de datos:
lmDiscrete =
DiscretizedBoundaryConditionData[SparseArray[, dof, 1],
SparseArray[, dof, dof], diriMat, diriRows,
diriVals, dof, 0, lmdof, 1];
El primer escaso array (que es cero) se agrega a la carga, el segundo (que también es cero en nuestro caso) se agrega a la matriz de rigidez. Luego se especifican diriMat, diriRows y diriVals. Por último, se dan algunos parámetros de tamaño. (En lo que tendré que pensar es en proporcionar una interfaz para generar DiscretizedBoundaryConditionData de tal manera que tales hacks no son correspondidos).
DeployBoundaryConditions[load, stiffness, lmDiscrete,
"ConstraintMethod" -> "Append"]
(*Solution,interpolation function and plot*)
solution = LinearSolve[stiffness, load];
(* Truncate solution back to dof size *)
solution = Take[solution, dof];
Cuando subes "MeshOptions" -> "MaxCellMeasure" -> 0.02, "MeshOrder" -> 2 obtienes algo como esto:
uif = ElementMeshInterpolation[mesh, solution];
Plot3D[uif[x1, x2], Element[x1, x2, reg], PlotRange -> All,
AxesLabel -> x1, x2]
Y
Plot[uif[0, y] - uif[1, y], y, 0, 1]
¿Todo está bien? 😉
La versión 11.0 tiene una condición de límite adicional PeriodicBoundaryCondition. Con esto puedes usar:
ufun = NDSolveValue[-Laplacian[u[x, y], x, y] == 10 x^5*y^3,
DirichletCondition[u[x, y] == 0, x == 0 && y == 0],
PeriodicBoundaryCondition[u[x, y], 0 <= y <= 1 && x == 1,
FindGeometricTransform[0, 0, 0, 1, 1, 0, 1, 1][[2]]],
PeriodicBoundaryCondition[u[x, y], 0 <= x <= 1 && y == 1,
FindGeometricTransform[0, 0, 1, 0, 0, 1, 1, 1][[
2]]], u, x, y [Element] Rectangle[0, 0, 1, 1]];
Plot3D[ufun[x, y], x, y [Element] Rectangle[0, 0, 1, 1]]
Esto también funciona en regiones no rectangulares. Una cosa a tener en cuenta es que solo necesitas una DirichletCondition ya que eso se asigna al lado opuesto y las dos condiciones ahora se asignan al otro lado opuesto, lo que las convierte en cuatro en total.
Reseñas y valoraciones del tutorial
Si entiendes que te ha resultado provechoso nuestro post, sería de mucha ayuda si lo compartes con más entusiastas de la programación de este modo contrubuyes a difundir nuestra información.