Skip to content

Commit 5a520eb

Browse files
yueminglyuemingl
committed
write to tecplot, fix integration weight, dirichlet boundary
1 parent 2f3be9b commit 5a520eb

File tree

9 files changed

+2735
-55
lines changed

9 files changed

+2735
-55
lines changed

‎double_hex3.1.dat‎

Lines changed: 1501 additions & 0 deletions
Large diffs are not rendered by default.

‎patch_triangle.grd‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,33 @@
1+
# UCD geometry file from Gridgen 15.08R1, a product of Pointwise, Inc.
2+
# 01-Nov-10 23:08:59
3+
#
4+
13 16 0 0 0
5+
1 0.000000000e+000 0.000000000e+000 0.000000000e+000
6+
2 1.000000000e+001 0.000000000e+000 0.000000000e+000
7+
3 1.000000000e+001 1.000000000e+001 0.000000000e+000
8+
4 0.000000000e+000 1.000000000e+001 0.000000000e+000
9+
5 5.000000000e+000 0.000000000e+000 0.000000000e+000
10+
6 1.000000000e+001 5.000000000e+000 0.000000000e+000
11+
7 5.000000000e+000 1.000000000e+001 0.000000000e+000
12+
8 0.000000000e+000 5.000000000e+000 0.000000000e+000
13+
9 5.000000000e+000 5.000000000e+000 0.000000000e+000
14+
10 7.500000000e+000 2.500000000e+000 0.000000000e+000
15+
11 2.500000000e+000 7.500000000e+000 0.000000000e+000
16+
12 7.500000000e+000 7.500000000e+000 0.000000000e+000
17+
13 2.500000000e+000 2.500000000e+000 0.000000000e+000
18+
1 0 tri 9 5 10
19+
2 0 tri 11 4 8
20+
3 0 tri 10 2 6
21+
4 0 tri 3 7 12
22+
5 0 tri 1 5 13
23+
6 0 tri 9 7 11
24+
7 0 tri 9 8 13
25+
8 0 tri 9 6 12
26+
9 0 tri 10 5 2
27+
10 0 tri 9 10 6
28+
11 0 tri 13 5 9
29+
12 0 tri 1 13 8
30+
13 0 tri 12 7 9
31+
14 0 tri 3 12 6
32+
15 0 tri 11 7 4
33+
16 0 tri 9 11 8

‎src/symjava/examples/Example6.java‎

Lines changed: 72 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,9 @@
44
import static symjava.math.SymMath.grad;
55
import static symjava.symbolic.Symbol.*;
66

7+
import java.util.HashMap;
78
import java.util.List;
9+
import java.util.Map;
810

911
import Jama.Matrix;
1012
import symjava.math.Transformation;
@@ -17,22 +19,41 @@
1719
import symjava.symbolic.Int;
1820
import symjava.symbolic.SymConst;
1921

22+
/**
23+
* Finite Element solver
24+
*
25+
*/
2026
public class Example6 {
2127
public static void main(String[] args) {
22-
//Read the mesh
23-
Mesh2D mesh = new Mesh2D("mesh1", x,y);
24-
mesh.readTriangleMesh("double_hex3.1.node", "double_hex3.1.ele");
25-
26-
//Our PDE equation
2728
Func u = new Func("u", x, y);
2829
Func v = new Func("v", x, y);
29-
Eq pde = new Eq(0.5*dot(grad(u), grad(v)) + 0.1*u*v, (x*x+y*y)*v);
30-
31-
solve(pde, mesh);
30+
31+
// //Our PDE equation
32+
// Eq pde = new Eq(0.5*dot(grad(u), grad(v)) + 0.1*u*v, (x*x+y*y)*v);
33+
// //Read the mesh
34+
// Mesh2D mesh = new Mesh2D("mesh1", x, y);
35+
// mesh.readTriangleMesh("double_hex3.1.node", "double_hex3.1.ele");
36+
// solve(pde, mesh, null, "double_hex3.1.dat");
37+
38+
//Another PDE equation with Dirichlet condition
39+
Eq pde2 = new Eq(dot(grad(u), grad(v)), (-2*(x*x+y*y)+36)*v);
40+
Mesh2D mesh2 = new Mesh2D("mesh2", x, y);
41+
//mesh2.readGridGenMesh("patch_triangle.grd");
42+
mesh2.readGridGenMesh("triangle.grd");
43+
//Mark boundary nodes
44+
double eps = 0.01;
45+
for(Node n : mesh2.nodes) {
46+
//if(Math.abs(1-Math.abs(n.coords[0]))<eps || Math.abs(1-Math.abs(n.coords[1]))<eps)
47+
if(Math.abs(3-Math.abs(n.coords[0]))<eps || Math.abs(3-Math.abs(n.coords[1]))<eps)
48+
n.setType(1);
49+
}
50+
Map<Integer, Double> diri = new HashMap<Integer, Double>();
51+
diri.put(1, 0.0);
52+
solve(pde2, mesh2, diri, "triangle.dat");
3253
}
3354

34-
public static void solve(Eq pde, Mesh2D mesh) {
35-
UnitRightTriangle tri = new UnitRightTriangle("Tri", r, s);
55+
public static void solve(Eq pde, Mesh2D mesh, Map<Integer, Double> dirichlet, String output) {
56+
System.out.println(String.format("PDE Weak Form: %s = %s", pde.lhs, pde.rhs));
3657

3758
//Create coordinate transformation
3859
SymConst x1 = new SymConst("x1");
@@ -63,6 +84,7 @@ public static void solve(Eq pde, Mesh2D mesh) {
6384
Expr sx = -jacMat[1][0]/jac; //sx = -yr/jac
6485
Expr sy = jacMat[0][0]/jac; //sy = xr/jac
6586

87+
UnitRightTriangle tri = new UnitRightTriangle("Tri", r, s);
6688
Int lhsInt[][] = new Int[shapeFuns.length][shapeFuns.length];
6789
Int rhsInt[] = new Int[shapeFuns.length];
6890
for(int i=0; i<shapeFuns.length; i++) {
@@ -109,9 +131,11 @@ public static void solve(Eq pde, Mesh2D mesh) {
109131
}
110132

111133
//Assemble the system
134+
System.out.println("Start assemble the system...");
135+
long begin = System.currentTimeMillis();
112136
List<Domain> eles = mesh.getSubDomains();
113137
double[][] matA = new double[mesh.nodes.size()][mesh.nodes.size()];
114-
double[][] vecb = new double[1][mesh.nodes.size()];
138+
double[] vecb = new double[mesh.nodes.size()];
115139
for(Domain d : eles) {
116140
Element e = (Element)d;
117141
double[] nodeCoords = e.getNodeCoords();
@@ -121,15 +145,47 @@ public static void solve(Eq pde, Mesh2D mesh) {
121145
int idxJ = e.nodes.get(j).index-1;
122146
matA[idxI][idxJ] += lhsNInt[i][j].eval(nodeCoords);
123147
}
124-
vecb[0][idxI] = rhsNInt[i].eval(nodeCoords);
148+
vecb[idxI] = rhsNInt[i].eval(nodeCoords);
125149
}
126150
}
151+
System.out.println("Assemble done! Time: "+(System.currentTimeMillis()-begin)+"ms");
127152

153+
System.out.println("Solving...");
154+
begin = System.currentTimeMillis();
128155
Matrix A = new Matrix(matA);
129-
Matrix b = new Matrix(vecb);
130-
Matrix x = A.solve(b);
131-
for(int i=0; i<x.getRowDimension(); i++) {
132-
System.out.println(x.get(i, 0));
156+
Matrix b = new Matrix(vecb, vecb.length);
157+
if(dirichlet != null) {
158+
for(Node n : mesh.nodes) {
159+
Double diri = dirichlet.get(n.type);
160+
if(diri != null) {
161+
setDirichlet(A, b, n.index-1, diri);
162+
}
163+
}
133164
}
165+
Matrix x = A.solve(b);
166+
// for(int i=0; i<x.getRowDimension(); i++) {
167+
// System.out.println(x.get(i, 0));
168+
// }
169+
mesh.writeTechplot(output, x.getArray());
170+
System.out.println("Solved! Time: "+(System.currentTimeMillis()-begin)+"ms");
171+
System.out.println("See the output file(Tecplot format) "+output+" for the solution.");
134172
}
173+
174+
public static void setDirichlet(Matrix A, Matrix b, int nodeIndex, double value) {
175+
int row = nodeIndex;
176+
int col = nodeIndex;
177+
A.set(row, col, 1.0);
178+
b.set(row, 0, value);
179+
for(int r=0; r<A.getRowDimension(); r++) {
180+
if(r != row) {
181+
A.set(r, col, 0.0);
182+
b.set(r, 0 , b.get(r, 0)-A.get(r, col)*value);
183+
}
184+
}
185+
for(int c=0; c<A.getColumnDimension(); c++) {
186+
if(c != col) {
187+
A.set(row, c, 0.0);
188+
}
189+
}
190+
}
135191
}

0 commit comments

Comments
 (0)