diff --git a/docs/_mesh/bd.md b/docs/_mesh/bd.md index 8a8dcb8..5507741 100644 --- a/docs/_mesh/bd.md +++ b/docs/_mesh/bd.md @@ -28,6 +28,8 @@ bdFlag = setboundary(node,elem,'Dirichlet','abs(x) + abs(y) == 1','Neumann','y== bdFlag = setboundary3(node,elem,'Dirichlet','(z==1) | (z==-1)'); ``` + + ## Local Labeling of Edges and Faces We label three edges of a triangle such that `bdFlag(t,i)` is the edge @@ -53,6 +55,76 @@ The ordering of edges is specified by `locEdge`. If we use `locEdge = [2 3; 1 3; +## Changing Ordering + +If we change the ordering of `elem`, the corresponding local faces are changed. Therefore when we sort the `elem`, we should sort the `bdFlag` simutaneously. For example, `sort(elem,2)` sorts the `elem` only, and leave `bdFlag` unchanged. Use `sortelem` (2D) or `sortelem3` (3D) to sort `elem` and `bdFlag` at the same time. + + +```matlab +[node,elem] = cubemesh([-1,1,-1,1,-1,1],2); +bdFlag = setboundary3(node,elem,'Dirichlet','x==1','Neumann','x~=1'); +figure(2); clf; +showmesh3(node,elem); +display(elem); display(bdFlag); +findnode3(node,[1,2,3,4,5,7,8]); +display('change to ascend ordering'); +[elem,bdFlag] = sortelem3(elem,bdFlag) +``` + + + elem = + + 1 2 3 7 + 1 4 3 7 + 1 5 6 7 + 1 5 8 7 + 1 2 6 7 + 1 4 8 7 + + + + + bdFlag = + 6x4 uint8 matrix + + 1 0 0 2 + 2 0 0 2 + 2 0 0 2 + 2 0 0 2 + 1 0 0 2 + 2 0 0 2 + + change to ascend ordering + + elem = + + 1 2 3 7 + 1 3 4 7 + 1 5 6 7 + 1 5 7 8 + 1 2 6 7 + 1 4 7 8 + + + + + bdFlag = + 6x4 uint8 matrix + + 1 0 0 2 + 2 0 0 2 + 2 0 0 2 + 2 0 2 0 + 1 0 0 2 + 2 0 2 0 + + +​ + +![png](mesh_figures/sc3doc_34_1.png) + + + ## Extract Boundary Edges and Faces We may extract boundary edges for a 2-D triangulation from `bdFlag` from the following code. If `elem` is sorted counterclockwise, the boundary edges inherits the orientation. See also `findboundary` and `findboundary3`. @@ -78,6 +150,80 @@ Or simply call [bdNode,bdEdge,isBdNode] = findboundary(elem,bdFlag); ``` +To find the outward normal direction of the boundary face, we use `gradbasis3` to get `Dlambda(t,:,k)` which is the gradient of $\lambda_k$. The outward normal direction of the kth face can be obtained by `-Dlambda(t,:,k)` which is independent of the ordering and orientation of `elem`. + + +```matlab +Dlambda = gradbasis3(node,elem); +T = auxstructure3(elem); +elem2face = T.elem2face; +face = T.face; +NF = size(face,1); +if ~isempty(bdFlag) + % Find Dirchelt boundary faces and nodes + isBdFace = false(NF,1); + isBdFace(elem2face(bdFlag(:,1) == 1,1)) = true; + isBdFace(elem2face(bdFlag(:,2) == 1,2)) = true; + isBdFace(elem2face(bdFlag(:,3) == 1,3)) = true; + isBdFace(elem2face(bdFlag(:,4) == 1,4)) = true; + DirichletFace = face(isBdFace,:); + % Find outwards normal direction of Neumann boundary faces + bdFaceOutDirec = zeros(NF,3); + bdFaceOutDirec(elem2face(bdFlag(:,1) == 2,1),:) = -Dlambda(bdFlag(:,1) == 2,:,1); + bdFaceOutDirec(elem2face(bdFlag(:,2) == 2,2),:) = -Dlambda(bdFlag(:,2) == 2,:,2); + bdFaceOutDirec(elem2face(bdFlag(:,3) == 2,3),:) = -Dlambda(bdFlag(:,3) == 2,:,3); + bdFaceOutDirec(elem2face(bdFlag(:,4) == 2,4),:) = -Dlambda(bdFlag(:,4) == 2,:,4); +end +% normalize the boundary face outwards direction +vl = sqrt(dot(bdFaceOutDirec,bdFaceOutDirec,2)); +idx = (vl==0); +NeumannFace = face(~idx,:); +bdFaceOutDirec(idx,:) = []; +vl(idx) = []; +bdFaceOutDirec = bdFaceOutDirec./[vl vl vl]; +display(DirichletFace); +display(NeumannFace); +display(bdFaceOutDirec); +``` + + + DirichletFace = + + 2x3 uint32 matrix + + 2 3 7 + 2 6 7 + + NeumannFace = + 10x3 uint32 matrix + + 1 2 3 + 1 2 6 + 1 3 4 + 1 4 8 + 1 5 6 + 1 5 8 + 3 4 7 + 4 7 8 + 5 6 7 + 5 7 8 + + + + + bdFaceOutDirec = + + 0 0 -1 + 0 -1 0 + 0 0 -1 + -1 0 0 + 0 -1 0 + -1 0 0 + 0 1 0 + 0 1 0 + 0 0 1 + 0 0 1 + @@ -112,10 +258,10 @@ display(bdFlag) ​ ​ bdFlag = ​ - 1 0 1 - 1 0 0 - 1 0 0 - 1 1 0 +​ 1 0 1 +​ 1 0 0 +​ 1 0 0 +​ 1 1 0 @@ -170,9 +316,9 @@ display(bdFlag) ​ ​ bdFlag = ​ - 0 1 0 0 - 0 0 0 0 - 1 0 0 0 +​ 0 1 0 0 +​ 0 0 0 0 +​ 1 0 0 0 @@ -188,3 +334,4 @@ flux). ## Remark It would save storage if we record boundary edges or faces only. The current data structure is convenient for the local refinement and coarsening since the boundary can be easily updated along with the change of elements. The matrix `bdFlag` is sparse but a dense matrix is used. We do not save `bdFlag` as a sparse matrix since updating sparse matrices is time-consuming. Instead we set up the type of `bdFlag` to `uint8` to minimize the waste of memory. + diff --git a/docs/_mesh/mesh_figures/sc3doc_31_1.png b/docs/_mesh/mesh_figures/sc3doc_31_1.png index 52b40f9..17cf7d8 100644 Binary files a/docs/_mesh/mesh_figures/sc3doc_31_1.png and b/docs/_mesh/mesh_figures/sc3doc_31_1.png differ diff --git a/docs/_mesh/sc.md b/docs/_mesh/sc.md index 554cf53..f05c780 100644 --- a/docs/_mesh/sc.md +++ b/docs/_mesh/sc.md @@ -38,8 +38,7 @@ findnode(node); ![png](mesh_figures/scdoc_3_0.png) -The basic data structure of a mesh consists of `node` and `elem`. -The integers `N, NE, NT, NE` represents the muber of vertice, edges, triangles, edges +The integers `N, NE, NT` represents the muber of vertice, edges, triangles respectively. The corresponding simplicial complex consists of vertices, edges and @@ -62,8 +61,8 @@ For indexing, ordering and orientation, there are always two versions: local and ## Indexing The indexing refers to the numbering of simplexes, e.g., which edge is -numbered as the first one. Each simplex in the simplicial complex has a unique index which is called the global index. In one triangle, the three vertices and three -edges have their local index from 1:3. +numbered as the first one. Each simplex in the simplicial complex has a unique index which is called the global index. Inside one triangle, the three vertices and three +edges have their local index `1:3`. In the assembling procedure of finite element methods, an element-wise matrix using the local index is firstly computed and then assembled to get a @@ -77,7 +76,7 @@ is sufficient. ### Pointers of vertices -The `NT x 3` matrix `elem` is indeed the pointer of the vertex indices. For example `elem(t,1)=25` means the first vertex of the triangle t is the 25-th vertex. +The `NT x 3` matrix `elem` is indeed the pointer of the vertex indices. For example `elem(t,1)=25` means the first vertex of `t` is the 25-th vertex. Similiary, the `NE x 2` matrix `edge` records the vertices pointer of edges. @@ -196,7 +195,7 @@ display('After rotation'); display(elem); display(bdFlag); ​ Ascend ordering will benefit the construction of local bases for high -order Lagrange elements or in general bases with orientation dependent. +order Lagrange elements and in general bases with orientation dependent. We may switch the default positive ordering of `elem` to the ascend ordering when generating data structure for finite element basis. However such @@ -263,7 +262,7 @@ the global ordering-orientation, it is the sign of the signed area In the output of `gradbasis`, `area` is always positive and an additional array `elemSign` is used to record the sign of the signed area. -`Dlambda(t,:,k)` is the gradient of the barycentric coordinate $\lambda_k$. Therefore the outward normal direction of the `k`th edge can be obtained by `-Dlambda(t,:,k)` which is independent of the ordering and orientation of triangle `t`. +`Dlambda(t,:,k)` is the gradient of the barycentric coordinate $\lambda_k$ associated to vertex $k$. Therefore the outward normal direction of the `k`th edge can be obtained by `-Dlambda(t,:,k)` which is independent of the orientation of triangle `t`. ### edge @@ -415,7 +414,7 @@ display(elem2edgeSign); 1 1 -1 - + ![png](mesh_figures/scdoc_16_1.png) @@ -460,6 +459,6 @@ display(elem2edgeSign); 1 -1 1 - + ![png](mesh_figures/scdoc_17_1.png) diff --git a/docs/_mesh/sc3.md b/docs/_mesh/sc3.md index aa75b73..77d957e 100644 --- a/docs/_mesh/sc3.md +++ b/docs/_mesh/sc3.md @@ -5,128 +5,112 @@ sidebar: nav: mesh --- -We describe the data structure of the simplicial complex associated to a three-dimensional triangulation give by `node,elem`. The `node` records -the coordinates of vertices and `elem` is the pointer from local to -global indices of vertices. See [Basic mesh data structure]({{ site.baseurl }}{% link _mesh/basic.md %}). +We describe the data structure of the simplicial complex associated to a three-dimensional triangulation give by `node,elem`. The `node` records the coordinates of vertices and `elem` is the pointer from local to global indices of vertices. See [Basic mesh data structure]({{ site.baseurl }}{% link _mesh/basic.md %}). A brief summary. - `edge`: ascending ordering, i.e. `edge(:,1) The multigrid solvers use the original ordering of `elem` obtained from either uniform refinement or bisection methods. So let `elemold=elem` before sort. +> The multigrid solvers use the original ordering of `elem` obtained from either uniform refinement or bisection methods. So let `elemold=elem` and use `elemold` in multigrid solvers. - Examples on the usage: `Poisson3RT0; Maxwell; Maxwell2;` + + ## Outline -The basic data structure of a mesh consists of `node` and `elem`. The corresponding simplicial complex consists of vertices, edges, faces, and tetrahedron. We shall discuss three issues +We describe the data structure of the simplicial complex associated to a +two-dimensional triangulation give by `node,elem`. The `node` records +the coordinates of vertices and `elem` is the pointer from the local to +the global indices of vertices. See [Basic mesh data structure]({{ site.baseurl }}{% link _mesh/basic.md %} ). + + +```matlab +node = [-1,-1,-1; 1,-1,-1; 1,1,-1; -1,1,-1; -1,-1,1; 1,-1,1; 1,1,1; -1,1,1]; +elem = [1,2,3,7; 1,6,2,7; 1,5,6,7; 1,8,5,7; 1,4,8,7; 1,3,4,7]; +clf; showmesh3(node,elem,[],'FaceAlpha',0.25); +view([-53,8]); +axis on +findelem3(node,elem); +findnode3(node); +``` + +![png](mesh_figures/meshbasicdoc_6_0.png) + +The integers `N, NE, NF, NT` represents the muber of vertice, edges, faces, tetrahedron +respectively. + +The corresponding simplicial complex consists of vertices, edges, faces, and +tetrahedron. We shall discuss the following three issues: - *Indexing* of simplexes -- *Ordering* of vertices -- *Orientation* of simplexes +- *Ordering* of vertices of simplexes +- *Orientatoin* of simplexes The indexing and ordering are related and the ordering and orientation -are mixed together. However, the indexing has nothing to do with the -orientation. The indexing and ordering are the combinatory structure, +are mixed together. However the indexing has nothing to do with the +orientation. The indexing and ordering are the combinarotry structure, i.e. only `elem` is needed, while the orientation also depends on `node`, -the geometry embedding of vertices. +the geometry emembdding of vertices. + +For indexing, ordering and orientation, there are always two versions: local and global. The pointer from the local and the global version is the most complicated issue. + -For indexing, ordering and orientation, there are always local and global versions. The relation between the local and global version is the most complicated issue. -## Indexing of Simplexes +## Indexing The indexing refers to the numbering of simplexes, e.g., which face is -numbered as the first one. There are two types of the indexing: local and -global. Each simplex in the simplicial complex has a unique index which -is called the global index. In one tetrahedron, the four vertices and four -faces have their local index from 1:4. +numbered as the first one. Each simplex in the simplicial complex has a unique index which is called the global index. Inside one tetrahedron, the four vertices and four +faces have their local index `1:4`. In the assembling procedure of finite element methods, an element-wise -matrix using the local indexing is first computed and then assembled to get a -big matrix using the global indexing. Thus the pointer from the local -indexing to the global indexing is indispensable. For bases independent of -the ordering and orientation, e.g., `P1` and `P2` elements, this pointer -is sufficient, otherwise, the inconsistency of the local ordering/orientation -and the global ordering/orientation should be taken into account. +matrix using the local index is firstly computed and then assembled to get a +big matrix using the global index. Thus the pointer from the local +index to the global one is indispensible. For bases independent of +the ordering and orientation, e.g., `P1` and `P2` elements, the pointer of indices +is sufficient. -### Local indexing +### Pointers of vertices -The tetrahedron consists of four vertices indexed as [1 2 3 4]. Each -tetrahedron contains four faces and six edges. They can be indexed as +The `NT x 4` matrix `elem` is indeed the pointer of the vertex indices. For example `elem(t,1)=25` means the first vertex of `t` is the 25-th vertex. - locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3]; - locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]; +Similiary, the `NE x 2` matrix `edge` records the vertices pointer of edges. + +### Local indexing of edges and faces + +The tetrahedron consists of four vertices indexed as `[1 2 3 4]`. Each +tetrahedron contains four faces and six edges. + +- Opposite indexing `locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3];` -In `locFace`, the $i$-th face is opposite to the $i$-th vertices and thus -this is called _opposite indexing_. In `locEdge`, it is the -_lexicographic indexing_ which is induced from the lexicographic ordering -of the six edges. The ordering of vertices of each face or edge will not +- Lexicographic indexing `locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4];` + +In `locFace`, the i-th face is opposite to the i-th vertices and thus called *opposite indexing*. In `locEdge`, it is the _lexicographic indexing_ which is induced from the lexicographic ordering of the six edges. The ordering of vertices inside each face or edge will not change the indexing. For example, the following `locFacec` and `locEdged` -has the same indexing as `locFace` and `locEdge` but a different ordering +has the same indexing as `locFace` and `locEdge` but with a different ordering of vertices. locFacec = [2 3 4; 1 4 3; 1 2 4; 1 3 2]; locEdge = [2 1; 3 1; 4 1; 3 2; 4 2; 4 3]; -Indeed any permutation of each simplex will represent the same simplex and -will not change the indexing. The ordering of vertices will affect the +Indeed any permutation of vertex indices will represent the same simplex and +will not change the indexing of the simplex. The ordering of vertices will affect the orientation and will be discussed later. -For a face consists of three vertices [1 2 3], there are two indexing +For a face consists of three vertices `[1 2 3]`, there are two indexing schemes of its three edges. -- Oppoiste indexing: `[2 3; 3 1; 1 2]` -- Lexicographic indexing: `[1 2; 1 3; 2 3]` +- Oppoiste indexing `[2 3; 3 1; 1 2]` +- Lexicographic indexing `[1 2; 1 3; 2 3]` Each indexing scheme has its advantage and disadvantage and which one -to chose depends on the consideration of ordering and orientation. - -### Global indexing and vertex pointers - -Each simplex in the simplicial complex has a unique index. It is -represented by vertices pointer from the local index to the global index -of vertices. - -The matrix `elem` is the pointer from local to global indices of vertices -of tetrahedron, e.g. `elem(t,1)=25` means the first vertex of the -tetrahedron t is the 25-th vertex. - -Similarly, the `NE x 2` matrix `edge` records all edges and the `NF x 3` by 3 -matrix `face` records all faces of the triangulation. These are vertices -pointers. We shall discuss the elementwise pointer from the local indices to -the global indices for edges and faces. - - - -```matlab -[node,elem] = cubemesh([-1,1,-1,1,-1,1],2); -showmesh3(node,elem,[],'FaceAlpha',0.25); -findelem3(node,elem); -findnode3(node,elem(:)); -display(elem); -``` +to use depends on the consideration of ordering and orientation. - - elem = - - 1 2 3 7 - 1 4 3 7 - 1 5 6 7 - 1 5 8 7 - 1 2 6 7 - 1 4 8 7 - - -![png](mesh_figures/sc3doc_6_1.png) - - -### Generate index pointers for edges and faces +### Global indexing of edges and faces One can easily collect edges and faces elementwise. The issue is the duplication. For example, each interior face will be counted twice. The @@ -176,10 +160,11 @@ display(face); 7 8 -​ - face = +​ - 18x3 uint32 matrix + face = + + 18x3 uint32 matrix 1 2 3 1 2 6 @@ -206,7 +191,7 @@ In iFEM, `N,NE,NF,NT` represents the number of vertices, edges, faces and tetrah N = size(node,1); NT = size(elem,1); NF = size(face,1); NE = size(edge,1); -In the assembling procedure, the matrix is always computed elementwise and then assemble to a big one. A pointer from the local index of a simplex to its global index is thus indispensable. +In the assembling procedure, the matrix is always computed elementwise and then assemble to a big one. A pointer from the local index to its global one is thus indispensable. **Elementwise pointers** @@ -214,7 +199,7 @@ In the assembling procedure, the matrix is always computed elementwise and then - `elem2face(1:NT, 1:4)` - `elem2edge(1:NT, 1:6)` -Such information is exactly stored in the output of `unique` function. For example, `elem2face(t,1) = 17` means the first face of t (spanned by `[2 3 4]`) is the 17-th element in the `face` matrix. +Such information is exactly stored in the output of `unique` function. For example, `elem2face(t,1) = 17` means the first face of `t` (spanned by `[2 3 4]`) is the `17`-th face in the `face` matrix. ```matlab @@ -226,22 +211,24 @@ display(elem2face); ``` - elem2edge = - - 6x6 uint32 matrix - - 1 2 6 8 10 12 - 3 2 6 11 13 12 - 4 5 6 15 16 18 - 4 7 6 17 16 19 - 1 5 6 9 10 18 - 3 7 6 14 13 19 +``` +elem2edge = + + 6x6 uint32 matrix + + 1 2 6 8 10 12 + 3 2 6 11 13 12 + 4 5 6 15 16 18 + 4 7 6 17 16 19 + 1 5 6 9 10 18 + 3 7 6 14 13 19 +``` + -​ elem2face = - 6x4 uint32 matrix + 6x4 uint32 matrix 13 5 3 1 15 5 6 4 @@ -254,15 +241,11 @@ display(elem2face); **Face to edge Pointer** -`face2edge(1:NF,1:3)` records the global indices of three edges of a -face. This pointer depends on the ordering of vertices of faces and the -indexing of local edges in a face. We list the following two important -cases. Other combinations are possible but not attractive. +`face2edge(1:NF,1:3)` records the global indices of three edges of a face. This pointer depends on the ordering of vertices of faces and the indexing of local edges. We list the following two important cases. Other combinations are possible but not attractive. - Ascend ordering. -All local faces and local edges are ascending ordered. - +All local faces and local edges are ascending ordered. locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3]; locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]; edgeOfFace = [1 2; 1 3; 2 3]; @@ -271,14 +254,14 @@ All local faces and local edges are ascending ordered. - Consistent ordering The local face is ordered such that the corresponding orientation is -consistent with the induced orientation. +consistent with the induced orientation from the tetrahedron. locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2]; locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]; edgeOfFace = [2 3; 3 1; 1 2]; locFace2edge = [6 5 4; 6 2 3; 5 3 1; 4 1 2]; -The global one can be obtained from the composition of `elem2face` and +The global pointer can be obtained from the composition of `elem2face` and `locFace2edge`. For example, for the ascending ordering scheme, face2edge(elem2face(:,1),:) = elem2edge(:,[4 5 6]); @@ -286,11 +269,11 @@ The global one can be obtained from the composition of `elem2face` and face2edge(elem2face(:,3),:) = elem2edge(:,[1 3 5]); face2edge(elem2face(:,4),:) = elem2edge(:,[1 2 4]); -## Ordering of Vertices -We discuss the ordering of vertices of simplexes. Again there are local -ordering and global ordering. They may not be consistent, and a sign array -is used to record the inconsistency if any. + +## Ordering + +We discuss the ordering of vertices of simplexes. The local and the global ordering may not be consistent and a sign array can be used to record the inconsistency. The local ordering refers to the ordering of vertices in `locFace` or `locEdge`, i.e. the ordering of the local index of vertices. For elements @@ -298,38 +281,29 @@ associated to faces or edges, the local ordering could be used in the formulation of the local basis and thus the ordering does matter. The global ordering refers to the ordering of vertices in `face` or -`edge`, i.e., the ordering of the global index of vertices. Note that -that in either local or global ordering, permutation of vertices will -represent the same simplex. To fix an ordering we need extra information. +`edge`, i.e., the ordering of the global index of vertices. Note that in either local or global ordering, permutation of vertices will represent the same simplex. To fix an ordering we need extra information. ### elem -The local ordering is always [1 2 3 4]. Any permutation of four vertices of a tetrahedron still represents the same tetrahedron. Such freedom provides a room to record more information like: +The local ordering is always `[1 2 3 4]`. Any permutation of four vertices of a tetrahedron still represents the same tetrahedron. Such freedom provides a room to record more information like: a global ordering of vertices, an orientation of element, and refinement rules (uniform refinement or bisection) etc. -* global ordering of vertices -* an orientation of element -* refinement rules (uniform refinement or bisection) +For 2-D triangulations, three vertices of a triangle in 2-D is sorted counter-cloclwise and the first vertex is chosen as the newest vertex. Such ordering enables the efficient implementation of local refinement and coarsening in 2-D; see [Bisection in Two Dimensions]({{ site.baseurl }}{% link _afem/bisect.md %}) and [Coarsening in Two Dimensions]({{ site.baseurl }}{% link _afem/coarsen.md %}). -For 2-D triangulations, three vertices of a triangle in 2-D is sorted counter-cloclwise and the first vertex is chosen as the newest vertex. Such ordering enables the efficient implementation of local refinement -and coarsening in 2-D; see [Bisection in Two Dimensions]({{ site.baseurl }}{% link _afem/bisect.md %}) -and [Coarsening in Two Dimensions]({{ site.baseurl }}{% link _afem/coarsen.md %}). +In 3-D, for the longest edge bisection, the newest vertex (with the maximal generation) is stored as the last (4-th) vertex of a tetrahedron. For [3-D Red Refinement]({{ site.baseurl }}{% link _mesh/uniformrefine3.md %}).), the ordering determines the shape regularity of refined triangulation. Permutation of vertices in `elem` could deteriorate the angle condition after the refinement. -In 3-D, for the longest edge bisection, the newest vertex (with the highest generation) is stored as the last (4-th) vertex of a tetrahedron. For [3-D Red Refinement]({{ site.baseurl }}{% link _mesh/uniformrefine3.md %}).), the ordering determines the shape regularity of refined triangulation. Permutation of vertices in `elem` could deteriorate the angle condition after the refinement. - -We shall reserve the ordering of `elem` from the mesh refinement and -coarsening since they are more subtle. We switch the ordering when -generating data structure for finite element basis and assemble the -matrix equation. Such sorting is hidden in the subroutines when a finite -element basis requiring ordering is generated. +We shall preserve the ordering of `elem` from the mesh refinement and +coarsening since they are more subtle. We switch the ordering inside the subroutines when generating data structure for finite element basis and assembleing the +matrix equation. Two types of ordering of `elem` is of particular importantance -- Ascend ordering -- Positive ordering +- ascend ordering +- positive ordering In the ascending ordering, the vertices of `elem` is sorted such that - elem(t,1) < elem(t,2) < elem(t,3) < elem(t,4). + elem(t,1) < elem(t,2) < elem(t,3) < elem(t,4) -Such ordering will benefit the construction of local bases for high order basis or basis with orientation. This can be easily achieved by `elem = sort(elem,2)`. One has to rotate the boundary flag accordingly. +Ascend ordering will benefit the construction of local bases for high +order Lagrange elements and in general bases with orientation dependent. This can be easily achieved by `elem = sort(elem,2)`. But to include correct boundary conditions, one has to rotate the boundary flag accordingly using `sortelem3`. ```matlab @@ -351,9 +325,7 @@ display(elem); In the positive ordering, the four vertices are ordered such that the -signed volume, the mix product of vectors (v12,v13,v14), is positive. -This is the default ordering used so far. `fixorder3` will switch the -vertices for elements with negative volume. +signed volume is positive. `fixorder3` will switch the vertices for elements with negative volume. ```matlab @@ -374,32 +346,29 @@ elem = fixorder3(node,elem) % switchs the vertices for elements with negative ### edge -For 3-D triangulations, we chose the ascending ordering both locally and globally. Namely +For 3-D triangulations, we chose the ascending ordering both locally and globally, i.e., locEdge(:,1) < locEdge(:,2); edge(:,1) < edge(:,2); -Recall that for `locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]`, it is ascend -ordered. The `edge` produced by `unique` function is also ascend -ordered. - -There might be inconsistency between the local and global ordering. -That is `edge(elem2edge(t,1),1)` may not be smaller than `edge(elem2edge(t,1),2)`. It will be more clear from the discussion of the corresponding orientation. +Recall that for `locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]`, it is ascending ordered. The `edge` produced by `unique` function is also ascending ordered. -For 2-D triangulations, the global ordering is still ascend ordered. But locally it may not. For example, for the consisitent ordering `locEdge = [2 3; 3 1; 1 2]`, then `locEdge(2,1) > locEdge(2,2)`. +There might be inconsistency between the local and global ordering. Namely `edge(elem2edge(t,1),1) > edge(elem2edge(t,1),2)` may happen. ### face For 3-D triangulations, the `face` produced by `unique` function is already sorted in the second dimension such that the global ordering is ascended i.e. `face(:,1) < face(:,2) < face(:,3)`. The local ordering in `locFace`, however, is not always ascend ordered. - locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3]; % Ascend ordering - locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2]; % Consistent ordering +```matlab +locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3]; % Ascend ordering +locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2]; % Consistent ordering +``` Again the local and the global ordering may not be consistent. That is face(elem2face(t,:),1) < face(elem2face(t,:),2) < face(elem2face(t,:),3) -may not be always true unless we use the ascending ordering in both `face` and `locFace`. +may not be always true unless the ascending ordering is used in both `face` and `locFace`. ## Orientation @@ -452,7 +421,7 @@ In the output of `gradbasis3`, `volume` is always positive and an additional array `elemSign` is used to record the sign of the signed volume. -`Dlambda(t,:,k)` is the gradient of $\lambda_k$ associated to vertex $k$. Therefore the outward normal direction of the $k$-th face is `-Dlambda(t,:,k)` which is independent of the ordering and orientation of the element. +`Dlambda(t,:,k)` is the gradient of $\lambda_k$ associated to vertex $k$. Therefore the outward normal direction of the k-th face is `-Dlambda(t,:,k)` which is independent of the orientation of element `t`. ### face @@ -488,18 +457,18 @@ elem2faceSign = reshape(sum(sign(diff(totalFace(:,[1:3,1]),1,2)),2),NT,4) -When both `elem` and `locFace` are ascend ordered, the orientation of the +When both `elem` and `locFace` are ascending ordered, the orientation of the global ordering is consistent with that of the local ordering. Thus -`elem2faceSign` is not needed for the ascending ordering in assembling the mass matrix. +`elem2faceSign` is not needed for the ascending ordering when assembling the mass matrix. -But for the ascending ordering system, an `elem2faceSign` will be used when assembling differential operators because the orientation for Stokes theorem is induced orientation. For example, when computing `div` operators on a positive orientated tetrahedron, the faces should be orientated by the outwards normal direction but the global faces may not be. +But for the ascending ordering system, an `elem2faceSign` will be used when assembling the differential operator because the orientation in the Stokes theorem is the induced orientation. For example, when computing `div` operators on a positive orientated tetrahedron, the faces should be orientated by the outwards normal direction but the global orientation of a face may not be. If `elem` is positive ordered and `locFace` is consistently ordered, then this inconsistency is already recorded in `elem2faceSign`. -For the ascending ordering of `elem` and `locFace`, we use $+1$ if the -orientation of a face is the same with the induced outwords normal direction in a certain elem, and $-1$ otherwise. Then the inconsistency is given by `elem2faceSign = [1 -1 1 -1]` by comparing +For the ascending ordering of `elem` and `locFace`, we use `1` if the +orientation of a face is the same with the induced outwords normal direction, and `-1` otherwise. Then the inconsistency is `elem2faceSign = [1 -1 1 -1]` by comparing -- The induced orientation: `locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2]`; +- the induced orientation: `locFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2]`; - the ascending orientation: `locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3]`. Here we use the *ascending orientation* to refer to the orientation given by the ascending ordering. @@ -513,8 +482,8 @@ In summary, The orientation of edges is simpler than faces. Globally we always chose the global ascend ordering orientation. Namely the orientation of an edge is from the vertex with the smaller index to the larger one. -Locally the local ascend ordering may not be consistent with the global one. See Lowest Order Edge Element. For the ascending ordering of `elem` and `locEdge`, the local and the global -orientation will be consistent and no `elem2edgeSign` is needed! +The global basis associated to this edge, however, depends only on its +global orientation. The array `elem2edgeSign(1:NT, 1:3)` records the inconsistency of a local ordering-orientation and the global orientation. ```matlab @@ -529,7 +498,7 @@ elem2edgeSign = reshape(direction,NT,6) elem2edgeSign = - 6�6 int8 matrix + 6x6 int8 matrix 1 1 1 1 1 1 1 1 1 1 1 1 @@ -538,13 +507,14 @@ elem2edgeSign = reshape(direction,NT,6) 1 1 1 -1 1 1 1 1 1 -1 1 1 - +For the ascending ordering of `elem` and `locEdge`, the local and the global +orientation will be consistent and no `elem2edgeSign` is needed! ### face to edge For the ascending ordering `edgeofFace = [1 2; 1 3; 2 3]`, the local and the global ordering is consistent and so is the ordering orientation. -Then it is not consistent with the induced positive (counter-clockwise) orientation of edges. When the edge direction is the same with the induced direction, we use sign $+1$, otherwise $-1$. Then `face2edgeSign = [+1 -1 +1]` records the inconsistency of the ascending orientation of the induced orientation. +Then it is not consistent with the induced positive (counter-clockwise) orientation of edges. When the edge direction is the same with the induced direction, we use sign `+1` otherwise `-1`. Then `face2edgeSign = [+1 -1 +1]` records the inconsistency of the ascending orientation of the induced orientation. For the consistent ordering, `edgeofFace = [2 3; 3 1; 1 2]` which is consistent with the induced positive orientation but then may not be consistent with the global orientation of edges. We construct `face2edgeSign` to record such inconsistency @@ -587,7 +557,7 @@ We summarize the two popular ordering and orientation schemes below. ### Ascend Odering and Orientation -The ascending ordering and orientation is more algebraic, determined by the indices of vertices. +The ascending ordering and orientation is more algebraic, determined by the indices of vertices not the coordinates. #### Ascend ordering The array `elem` is sorted such that @@ -607,9 +577,8 @@ Then due to the ascending ordering of `elem`, globally the `edge` and `face` als One can easily see the benefit: the ordering of local edges and local faces is consistent with the global ones and so is their corresponding orientation. -#### Orientation -We chose the global ordering orientation for each element. We chose the orientation corresponding to the ascending ordering for edges -and faces. That is +#### Ordering Orientation +We chose the global ordering orientation for each element. * `elem: sign(v12,v13,v14)` * `face`: the normal vector is given by `cross(v12,v13)` @@ -620,32 +589,37 @@ For faces and edges, the orientation of the ascending ordering and the induced o * `elem2faceSign = [1 -1 1 -1];` * `face2edgeSign = [1 -1 1];` + + ### Positive Ordering and Orientation -The positive orientation and ordering is more geometrically consistent in the sense that the orientation of an element is locally consistent with the orientation of the local boundary faces. But it introduces inconsistency with the global orientation of a simplex. +The positive orientation and ordering is geometrically consistent in the sense that the orientation of an element is locally consistent with the orientation of the local boundary faces. But it introduces inconsistency with the global orientation of a simplex. #### Positive and consistent ordering -The vertices of `elem` is sorted such that the signed volume is always positive, i.e. the four vertices follows the right hand rule. +The vertices of `elem` is sorted such that the signed volume is always positive, i.e., the four vertices follows the right hand rule. The four faces of a tetrahedron are ordered consistently as locFace = [2 3 4; 1 3 4; 1 2 4; 1 2 3]; -The six edges of a tetrahedron still ascend ordered +The six edges of a tetrahedron are still ascending ordered locEdge = [1 2; 1 3; 1 4; 2 3; 2 4; 3 4]; -Three edges of a face is ordered consistently +Three edges of a face is ordered consistently not ascending edgeofFace = [2 3; 3 1; 1 2]; #### Orientation -the ascending ordering orientation is used for the global orientation of `edge` and `face` arrays. The inconsistency of the local and the global orientation is recorded in +The ascending ordering orientation is used for the global orientation of `edge` and `face` arrays. The inconsistency of the local and the global orientation is recorded in `elem2faceSign` and `elem2edgeSign`. -### An Example + + +### Example + We show two tetrahedron with the ascending ordering. @@ -684,10 +658,9 @@ elem2face = uint32(reshape(jf,NT,4)) 1 4 5 7 -​ - edge = - - 9x2 uint32 matrix + + edge = + 9x2 uint32 matrix 1 4 1 5 @@ -700,10 +673,9 @@ elem2face = uint32(reshape(jf,NT,4)) 5 8 -​ + face = - - 7x3 uint32 matrix + 7x3 uint32 matrix 1 4 5 1 4 7 @@ -716,32 +688,28 @@ elem2face = uint32(reshape(jf,NT,4)) ​ elem2edge = - - 2x6 uint32 matrix + 2x6 uint32 matrix 1 2 4 5 7 9 1 2 3 5 6 8 -​ + elem2face = - 2x4 uint32 matrix 7 5 3 1 6 4 2 1 -​ + v = - 0.6667 0.6667 -​ + elemSign = - -1 1 @@ -757,152 +725,3 @@ elem2face = uint32(reshape(jf,NT,4)) Since we are using the ascending ordering, the inconsistency with the induced orientation is * `elem2faceSign = [1 -1 1 -1];` * `face2edgeSign = [1 -1 1];` - -## Boundary Faces and Boundary Conditions - -We use `bdFlag` to record the boundary condition; see [Data Structure: Boundary Conditions]({{ site.baseurl }}{% link _mesh/bd.md %}).) for details. For short, `bdFlag` has the same size with `elem`, and records the boundary type of each local faces. If we change the ordering of `elem`, the corresponding local faces are changed. Therefore when we sort the `elem`, we should sort the `bdFlag` respectively. We use `sortelem3` to sort `elem` and `bdFlag` at the same time. Note that `sort(elem,2)` sorts the `elem` only, and leave -`bdFlag` unchanged. - - -```matlab -[node,elem] = cubemesh([-1,1,-1,1,-1,1],2); -bdFlag = setboundary3(node,elem,'Dirichlet','x==1','Neumann','x~=1'); -figure(2); clf; -showmesh3(node,elem); -display(elem); display(bdFlag); -findnode3(node,[1,2,3,4,5,7,8]); -display('change to ascend ordering'); -[elem,bdFlag] = sortelem3(elem,bdFlag) -``` - - - elem = - - 1 2 3 7 - 1 4 3 7 - 1 5 6 7 - 1 5 8 7 - 1 2 6 7 - 1 4 8 7 - - -​ - bdFlag = - - 6x4 uint8 matrix - - 1 0 0 2 - 2 0 0 2 - 2 0 0 2 - 2 0 0 2 - 1 0 0 2 - 2 0 0 2 - - change to ascend ordering - - elem = - - 1 2 3 7 - 1 3 4 7 - 1 5 6 7 - 1 5 7 8 - 1 2 6 7 - 1 4 7 8 - - -​ - bdFlag = - - 6x4 uint8 matrix - - 1 0 0 2 - 2 0 0 2 - 2 0 0 2 - 2 0 2 0 - 1 0 0 2 - 2 0 2 0 - - - - -​ - -![png](mesh_figures/sc3doc_34_1.png) - - -We can use `bdFlag` to find the boundary nodes, edges and faces. To find the outward normal direction of the boundary face, we use `gradbasis3` to get `Dlambda(t,:,k)` which is the gradient of $\lambda_k$. The outward normal direction of the kth face can be obtained by `-Dlambda(t,:,k)` which is independent of the ordering and orientation of `elem`. - - -```matlab -Dlambda = gradbasis3(node,elem); -T = auxstructure3(elem); -elem2face = T.elem2face; -face = T.face; -NF = size(face,1); -if ~isempty(bdFlag) - % Find Dirchelt boundary faces and nodes - isBdFace = false(NF,1); - isBdFace(elem2face(bdFlag(:,1) == 1,1)) = true; - isBdFace(elem2face(bdFlag(:,2) == 1,2)) = true; - isBdFace(elem2face(bdFlag(:,3) == 1,3)) = true; - isBdFace(elem2face(bdFlag(:,4) == 1,4)) = true; - DirichletFace = face(isBdFace,:); - % Find outwards normal direction of Neumann boundary faces - bdFaceOutDirec = zeros(NF,3); - bdFaceOutDirec(elem2face(bdFlag(:,1) == 2,1),:) = -Dlambda(bdFlag(:,1) == 2,:,1); - bdFaceOutDirec(elem2face(bdFlag(:,2) == 2,2),:) = -Dlambda(bdFlag(:,2) == 2,:,2); - bdFaceOutDirec(elem2face(bdFlag(:,3) == 2,3),:) = -Dlambda(bdFlag(:,3) == 2,:,3); - bdFaceOutDirec(elem2face(bdFlag(:,4) == 2,4),:) = -Dlambda(bdFlag(:,4) == 2,:,4); -end -% normalize the boundary face outwards direction -vl = sqrt(dot(bdFaceOutDirec,bdFaceOutDirec,2)); -idx = (vl==0); -NeumannFace = face(~idx,:); -bdFaceOutDirec(idx,:) = []; -vl(idx) = []; -bdFaceOutDirec = bdFaceOutDirec./[vl vl vl]; -display(DirichletFace); -display(NeumannFace); -display(bdFaceOutDirec); -``` - - - DirichletFace = - - 2x3 uint32 matrix - - 2 3 7 - 2 6 7 - - -​ - NeumannFace = - - 10x3 uint32 matrix - - 1 2 3 - 1 2 6 - 1 3 4 - 1 4 8 - 1 5 6 - 1 5 8 - 3 4 7 - 4 7 8 - 5 6 7 - 5 7 8 - - -​ - bdFaceOutDirec = - - 0 0 -1 - 0 -1 0 - 0 0 -1 - -1 0 0 - 0 -1 0 - -1 0 0 - 0 1 0 - 0 1 0 - 0 0 1 - 0 0 1 - diff --git a/solver/amgMaxwellinterface.m b/solver/amgMaxwellinterface.m new file mode 100644 index 0000000..de63a84 --- /dev/null +++ b/solver/amgMaxwellinterface.m @@ -0,0 +1,256 @@ +function [x,info] = amgMaxwellinterface(A,b,node,edge,option) +%% AMGMAXWELL algebraic multigrid solver for Maxwell equations. +% +% x = amgMaxwell(A,b,node,edge) attempts to solve the system of +% linear equations Ax = b using multigrid type solver. The linear system is +% obtained by either the first or second family of linear edge element +% discretization of the Maxwell equation; See doc Maxwell. +% +% amgMaxwell is more algebraic than mgMaxwell but still requires geometric +% information node and edge. Grapha Laplacian of vertices are used as +% auxiliary Poisson operator and amg is used as Poisson solver. +% +% Input +% - A: curl(alpha curl) + beta I +% - b: right hand side +% - node,edge: mesh information +% - options: extra structures +% +% By default, the HX preconditioned PCG is used which works well for +% symmetric positive definite matrices (e.g. arising from eddy current +% simulation). For symmetric indefinite matrices (e.g. arising from time +% harmonic Maxwell equation), set option.solver = 'minres' or 'bicg' or +% 'gmres' to try other Krylov method with HX preconditioner. +% +% See also mg, mgMaxwell, mgMaxwellsaddle +% +% Reference: +% R. Hiptmair and J. Xu, Nodal Auxiliary Space Preconditioning in +% H(curl) and H(div) Spaces. SIAM J. Numer. Anal., 45(6):2483-2509, 2007. +% +% Copyright (C) Long Chen. See COPYRIGHT.txt for details. + +t = cputime; +%% Initial check +Ndof = length(b); % number of dof +N = size(node,1); % number of nodes +NE = size(edge,1); % number of edge; +dim = size(node,2); +% Assign default values to unspecified parameters +if ~exist('option','var') + option = []; +end +option = mgoptions(option,length(b)); % parameters +x0 = option.x0; +% N0 = option.N0; +tol = option.tol; +% mu = option.smoothingstep; preconditioner = option.preconditioner; coarsegridsolver = option.coarsegridsolver; +printlevel = option.printlevel; %setupflag = option.setupflag; +maxIt = 400; % increase the default step (200) for Maxwell equations +% tol = 1e-7; % reset the tol +% Check for zero right hand side +if (norm(b) == 0) % if rhs vector is all zeros + x = zeros(Ndof,1); % then solution is all zeros + flag = 0; + itStep = 0; + err = 0; + time = cputime - t; + info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:))); + return +end + +%% Transfer operators from nodal element to edge element +if isfield(option,'isBdEdge') + isBdEdge = option.isBdEdge; +else + deg = sum(spones(A(1:NE,1:NE)),2); + isBdEdge = (deg == 1); +end +if Ndof == NE % lowest order edge element + II = node2edgematrix(node,edge,isBdEdge); +elseif Ndof >= 2*NE % first or second order edge element + II = node2edgematrix1(node,edge,isBdEdge); +end +IIt = II'; +[grad,isBdNode] = gradmatrix(edge,isBdEdge); +gradt = grad'; + +%% Free node +% isFreeNode = false(N,1); +% isFreeNode(edge(isBdEdge,:)) = true; +% isFreeNode = true(N,1); + +%% Auxiliary Poisson matrix +% - A: curl(alpha curl) + beta I +% - AP: - div(alpha grad) + beta I +% - BP: - div(beta grad) + +if isfield(option,'AP') + AP = option.AP; +else + % build graph Laplacian to approximate AP + edgeVec = node(edge(:,2),:) - node(edge(:,1),:); + edgeLength = sqrt(sum(edgeVec.^2,2)); + if isfield(option,'alpha') % resacle by the magnetic coefficients + if isreal(option.alpha) && (length(option.alpha) == NE) + alpha = option.alpha; + else % option.alpha is a function + edgeMiddle = (node(edge(:,2),:) + node(edge(:,1),:))/2; + alpha = option.alpha(edgeMiddle); + end + else + alpha = 1; + end + % edge weight: h*alpha + % AP = gradt*spdiags(edgeLength.*alpha,0,NE,NE)*grad; + i1 = (1:NE)'; j1 = double(edge(:,1)); s1 = sqrt(edgeLength.*alpha); + i2 = (1:NE)'; j2 = double(edge(:,2)); s2 = -s1; + isFreeEdge = ~isBdEdge; + G = sparse([i1(isFreeEdge);i2(isFreeEdge)],... + [j1(isFreeEdge);j2(isFreeEdge)],... + [s1(isFreeEdge);s2(isFreeEdge)],NE,N); + AP = G'*G; + % lumped mass matrix: h^3*beta + if isfield(option,'beta') % resacle by the dielectric coefficients + if isreal(option.beta) && (length(option.beta) == NE) + beta = option.beta; + else % option.beta is a function + edgeMiddle = (node(edge(:,2),:) + node(edge(:,1),:))/2; + beta = option.beta(edgeMiddle); + end + else + beta = 1; + end + M = accumarray(edge(:),repmat((edgeLength.^3).*beta,2,1),[N 1]); + AP = AP + spdiags(M,0,N,N); +end +% BP is Galerkin projection to the free node space +% boundary nodes +bdidx = zeros(N,1); +if isempty(isBdNode) % Neumann boundary condition + bdidx = 1e-6; +else + bdidx(isBdNode) = 1; +end +Tbd = spdiags(bdidx,0,N,N); +BP = gradt*A(1:NE,1:NE)*grad + Tbd; +% BP = gradt*spdiags(edgeLength.*beta,0,NE,NE)*grad + Tbd; + +%% Transfer operators between multilevel meshes +setupOption.solver = 'NO'; +[x,info,APi,Ri,RRi,ResAP,ProAP] = amg(AP,ones(N,1),setupOption); %#ok +[x,info,BPi,Si,SSi,ResBP,ProBP] = amg(BP,ones(N,1),setupOption); %#ok +D = diag(A); + +%% Krylov iterative methods with HX preconditioner +k = 1; +err = 1; +switch upper(option.outsolver) + case 'CG' + if printlevel>=1 + fprintf('Conjugate Gradient Method using HX preconditioner \n'); + end + x = x0; + r = b - A*x; + nb = norm(b); + err = zeros(maxIt,2); + err(1,:) = norm(r)/nb; + while (max(err(k,:)) > tol) && (k <= maxIt) + % compute Br by HX preconditioner + Br = HXpreconditioner(r); + % update tau, beta, and p + rho = r'*Br; % r'*Br = e'*ABA*e approximates e'*A*e + if k==1 + p = Br; + else + beta = rho/rho_old; + p = Br + beta*p; + end + % update alpha, x, and r + Ap = A*p; + alpha = rho/(Ap'*p); + r = r - alpha*Ap; + x = x + alpha*p; + rho_old = rho; + % compute err for the stopping criterion + k = k + 1; + err(k,1) = sqrt(abs(rho/(x'*b))); % approximate relative error in energy norm + err(k,2) = norm(r)/nb; % relative error of the residual in L2-norm + if printlevel >= 2 + fprintf('#dof: %8.0u, HXCG iter: %2.0u, err = %12.8g\n',... + Ndof, k, max(err(k,:))); + end + end + err = err(1:k,:); + itStep = k-1; + if k > maxIt || (max(err(end,:))>tol) + flag = 1; + else + flag = 0; + end + case 'MINRES' + fprintf('Minimum Residual Method with HX preconditioner \n') + [x,flag,err,itStep] = minres(A,b,tol,maxIt,@HXpreconditioner,[],x0); +% x = minres(A,b,tol,maxIt,@HXpreconditioner,[],x0); + case 'GMRES' + fprintf('General Minimum Residual Method with HX preconditioner \n') + [x,flag,err,itStep] = gmres(A,b,10,tol,maxIt,@HXpreconditioner,[],x0); + itStep = (itStep(1) -1)*10 + itStep(2); +end + +%% Output +time = cputime - t; +if printlevel >= 1 + fprintf('#dof: %8.0u, #nnz: %8.0u, iter: %2.0u, err = %8.4e, time = %4.2g s\n',... + Ndof, nnz(A), itStep, max(err(end,:)), time) +end +if (flag == 1) && (printlevel>0) + fprintf('NOTE: the iterative method does not converge'); +end +info = struct('solverTime',time,'itStep',itStep,'error',err,'flag',flag,'stopErr',max(err(end,:))); + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% subfunctions HXpreconditioner +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + function Br = HXpreconditioner(r) + %% 1. Smoothing in the finest grid of the original system +% eh = triu(A)\(D.*(tril(A)\r)); % Gauss-Seidal. less iteration steps + eh = 0.75*r./D; % Jacobi method. less computational time + + %% 2. Correction in the auxiliary spaces + % Part1: II*(AP)^{-1}*II^t + rc = reshape(IIt*r(1:size(IIt,2)),N,dim); % transfer to the nodal linear element space + level = size(APi,1); + ri = cell(level,1); + ei = cell(level,1); + ri{level} = rc; + for i = level:-1:2 + ei{i} = Ri{i}\ri{i}; + ri{i-1} = ResAP{i}*(ri{i}-APi{i}*ei{i}); + end + ei{1} = APi{1}\ri{1}; + for i = 2:level + ei{i} = ei{i} + ProAP{i-1}*ei{i-1}; + ei{i} = ei{i} + RRi{i}\(ri{i} - APi{i}*ei{i}); + end + eaux = [II*reshape(ei{level},dim*N,1); zeros(Ndof-size(IIt,2),1)]; % transfer back to the edge element space + % Part2: grad*(BP)^{-1}*grad^t + level = size(BPi,1); + ri{level} = gradt*r(1:NE); % transfer the residual to the null space + for i = level:-1:2 + ei{i} = Si{i}\ri{i}; + ri{i-1} = ResBP{i}*(ri{i}-BPi{i}*ei{i}); + end + ei{1} = BPi{1}\ri{1}; + for i=2:level + ei{i} = ei{i} + ProBP{i-1}*ei{i-1}; + ei{i} = ei{i} + SSi{i}\(ri{i}-BPi{i}*ei{i}); + end + eaux = eaux + [grad*ei{level}; zeros(Ndof-NE,1)]; % transfer back to the edge element space + % combined correction in the finest grid + Br = eh + eaux; + end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +end