Skip to content

Commit

Permalink
update sc and sc3
Browse files Browse the repository at this point in the history
  • Loading branch information
@lyc102
lyc102 committed Aug 30, 2021
1 parent 912c396 commit f40de46
Show file tree
Hide file tree
Showing 5 changed files with 569 additions and 348 deletions.
161 changes: 154 additions & 7 deletions docs/_mesh/bd.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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`.
Expand All @@ -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




Expand Down Expand Up @@ -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



Expand Down Expand Up @@ -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



Expand All @@ -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.

Binary file modified docs/_mesh/mesh_figures/sc3doc_31_1.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
17 changes: 8 additions & 9 deletions docs/_mesh/sc.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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.

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -415,7 +414,7 @@ display(elem2edgeSign);
1 1 -1




![png](mesh_figures/scdoc_16_1.png)

Expand Down Expand Up @@ -460,6 +459,6 @@ display(elem2edgeSign);
1 -1 1




![png](mesh_figures/scdoc_17_1.png)
Loading

0 comments on commit f40de46

Please sign in to comment.