Skip to content

Commit a43b8e6

Browse files
rezgarshakerirezgarshakeri
committed
Tried with stretched element, got same convergence rate.
1 parent efd09de commit a43b8e6

34 files changed

Lines changed: 400 additions & 160 deletions

‎examples/Hdiv-mass/conv_test.sh‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,8 @@ declare -A run_flags
3333
run_flags[problem]=mass2d
3434
run_flags[dm_plex_dim]=$dim
3535
run_flags[dm_plex_box_faces]=2,2
36+
run_flags[dm_plex_box_lower]=0,0
37+
run_flags[dm_plex_box_upper]=1,0.1
3638
else
3739
run_flags[problem]=mass3d
3840
run_flags[dm_plex_dim]=$dim
@@ -42,7 +44,7 @@ declare -A run_flags
4244
declare -A test_flags
4345
test_flags[res_start]=2
4446
test_flags[res_stride]=1
45-
test_flags[res_end]=5
47+
test_flags[res_end]=10
4648

4749
file_name=conv_test_result.csv
4850

‎examples/Hdiv-mass/conv_test_result.csv‎

Lines changed: 9 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,10 @@
11
run,mesh_res,error_u
2-
0,2,0.55094
3-
1,3,0.23091
4-
2,4,0.12473
5-
3,5,0.07812
2+
0,2,0.06228
3+
1,3,0.02836
4+
2,4,0.01616
5+
3,5,0.01043
6+
4,6,0.00728
7+
5,7,0.00537
8+
6,8,0.00413
9+
7,9,0.00327
10+
8,10,0.00265

‎examples/Hdiv-mass/convrate_mass.png‎

-2.11 KB
Loading

‎examples/Hdiv-mass/qfunctions/poisson-error2d.h‎

Lines changed: 6 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@
2121
#define ERROR2D_H
2222

2323
#include <math.h>
24-
24+
#include "utils.h"
2525
// -----------------------------------------------------------------------------
2626
// Compuet error
2727
// -----------------------------------------------------------------------------
@@ -41,17 +41,13 @@ CEED_QFUNCTION(SetupError2D)(void *ctx, const CeedInt Q,
4141
// Setup, J = dx/dX
4242
const CeedScalar J[2][2] = {{dxdX[0][0][i], dxdX[1][0][i]},
4343
{dxdX[0][1][i], dxdX[1][1][i]}};
44-
const CeedScalar detJ = J[0][0]*J[1][1] - J[0][1]*J[1][0];
44+
const CeedScalar det_J = MatDet2x2(J);
4545
// Compute Piola map:uh = J*u/detJ
46-
CeedScalar uh[2];
47-
for (CeedInt k = 0; k < 2; k++) {
48-
uh[k] = 0;
49-
for (CeedInt m = 0; m < 2; m++)
50-
uh[k] += J[k][m] * u[m][i]/detJ;
51-
}
46+
CeedScalar u1[2] = {u[0][i], u[1][i]}, uh[2];
47+
AlphaMatVecMult2x2(1/det_J, J, u1, uh);
5248
// Error
53-
error[i+0*Q] = (uh[0] - target[i+0*Q])*(uh[0] - target[i+0*Q])*w[i]*detJ;
54-
error[i+1*Q] = (uh[1] - target[i+1*Q])*(uh[1] - target[i+1*Q])*w[i]*detJ;
49+
error[i+0*Q] = (uh[0] - target[i+0*Q])*(uh[0] - target[i+0*Q])*w[i]*det_J;
50+
error[i+1*Q] = (uh[1] - target[i+1*Q])*(uh[1] - target[i+1*Q])*w[i]*det_J;
5551

5652
} // End of Quadrature Point Loop
5753

‎examples/Hdiv-mass/qfunctions/poisson-error3d.h‎

Lines changed: 7 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -21,23 +21,7 @@
2121
#define ERROR3D_H
2222

2323
#include <math.h>
24-
25-
// -----------------------------------------------------------------------------
26-
// Compute determinant of 3x3 matrix
27-
// -----------------------------------------------------------------------------
28-
#ifndef DetMat
29-
#define DetMat
30-
CEED_QFUNCTION_HELPER CeedScalar ComputeDetMat(const CeedScalar A[3][3]) {
31-
// Compute det(A)
32-
const CeedScalar B11 = A[1][1]*A[2][2] - A[1][2]*A[2][1];
33-
const CeedScalar B12 = A[0][2]*A[2][1] - A[0][1]*A[2][2];
34-
const CeedScalar B13 = A[0][1]*A[1][2] - A[0][2]*A[1][1];
35-
CeedScalar detA = A[0][0]*B11 + A[1][0]*B12 + A[2][0]*B13;
36-
37-
return detA;
38-
};
39-
#endif
40-
24+
#include "utils.h"
4125
// -----------------------------------------------------------------------------
4226
// Compuet error
4327
// -----------------------------------------------------------------------------
@@ -58,18 +42,14 @@ CEED_QFUNCTION(SetupError3D)(void *ctx, const CeedInt Q,
5842
const CeedScalar J[3][3] = {{dxdX[0][0][i], dxdX[1][0][i], dxdX[2][0][i]},
5943
{dxdX[0][1][i], dxdX[1][1][i], dxdX[2][1][i]},
6044
{dxdX[0][2][i], dxdX[1][2][i], dxdX[2][2][i]}};
61-
const CeedScalar detJ = ComputeDetMat(J);
45+
const CeedScalar det_J = MatDet3x3(J);
6246
// Compute Piola map:uh = J*u/detJ
63-
CeedScalar uh[3];
64-
for (CeedInt k = 0; k < 3; k++) {
65-
uh[k] = 0;
66-
for (CeedInt m = 0; m < 3; m++)
67-
uh[k] += J[k][m] * u[m][i]/detJ;
68-
}
47+
CeedScalar u1[3] = {u[0][i], u[1][i], u[2][i]}, uh[3];
48+
AlphaMatVecMult3x3(1/det_J, J, u1, uh);
6949
// Error
70-
error[i+0*Q] = (uh[0] - target[i+0*Q])*(uh[0] - target[i+0*Q])*w[i]*detJ;
71-
error[i+1*Q] = (uh[1] - target[i+1*Q])*(uh[1] - target[i+1*Q])*w[i]*detJ;
72-
error[i+2*Q] = (uh[2] - target[i+2*Q])*(uh[2] - target[i+2*Q])*w[i]*detJ;
50+
error[i+0*Q] = (uh[0] - target[i+0*Q])*(uh[0] - target[i+0*Q])*w[i]*det_J;
51+
error[i+1*Q] = (uh[1] - target[i+1*Q])*(uh[1] - target[i+1*Q])*w[i]*det_J;
52+
error[i+2*Q] = (uh[2] - target[i+2*Q])*(uh[2] - target[i+2*Q])*w[i]*det_J;
7353

7454
} // End of Quadrature Point Loop
7555

‎examples/Hdiv-mass/qfunctions/poisson-mass2d.h‎

Lines changed: 8 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@
2121
#define POISSON_MASS2D_H
2222

2323
#include <math.h>
24-
24+
#include "utils.h"
2525
// -----------------------------------------------------------------------------
2626
// This QFunction applies the mass operator for a vector field of 2 components.
2727
//
@@ -54,25 +54,19 @@ CEED_QFUNCTION(SetupMass2D)(void *ctx, CeedInt Q, const CeedScalar *const *in,
5454
// Setup, J = dx/dX
5555
const CeedScalar J[2][2] = {{dxdX[0][0][i], dxdX[1][0][i]},
5656
{dxdX[0][1][i], dxdX[1][1][i]}};
57-
const CeedScalar detJ = J[0][0]*J[1][1] - J[0][1]*J[1][0];
57+
const CeedScalar det_J = MatDet2x2(J);
5858

59-
const CeedScalar u1[2] = {u[0][i], u[1][i]};
59+
CeedScalar u1[2] = {u[0][i], u[1][i]}, v1[2];
6060
// *INDENT-ON*
6161
// Piola map: J^T*J*u*w/detJ
6262
// 1) Compute J^T * J
63-
CeedScalar JTJ[2][2];
64-
for (CeedInt j = 0; j < 2; j++) {
65-
for (CeedInt k = 0; k < 2; k++) {
66-
JTJ[j][k] = 0;
67-
for (CeedInt m = 0; m < 2; m++)
68-
JTJ[j][k] += J[m][j] * J[m][k];
69-
}
70-
}
63+
CeedScalar JT_J[2][2];
64+
AlphaMatTransposeMatMult2x2(1., J, J, JT_J);
65+
7166
// 2) Compute J^T*J*u * w /detJ
67+
AlphaMatVecMult2x2(w[i]/det_J, JT_J, u1, v1);
7268
for (CeedInt k = 0; k < 2; k++) {
73-
v[k][i] = 0;
74-
for (CeedInt m = 0; m < 2; m++)
75-
v[k][i] += JTJ[k][m] * u1[m] * w[i]/detJ;
69+
v[k][i] = v1[k];
7670
}
7771
} // End of Quadrature Point Loop
7872

‎examples/Hdiv-mass/qfunctions/poisson-mass3d.h‎

Lines changed: 7 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -21,23 +21,7 @@
2121
#define POISSON_MASS3D_H
2222

2323
#include <math.h>
24-
25-
// -----------------------------------------------------------------------------
26-
// Compute determinant of 3x3 matrix
27-
// -----------------------------------------------------------------------------
28-
#ifndef DetMat
29-
#define DetMat
30-
CEED_QFUNCTION_HELPER CeedScalar ComputeDetMat(const CeedScalar A[3][3]) {
31-
// Compute det(A)
32-
const CeedScalar B11 = A[1][1]*A[2][2] - A[1][2]*A[2][1];
33-
const CeedScalar B12 = A[0][2]*A[2][1] - A[0][1]*A[2][2];
34-
const CeedScalar B13 = A[0][1]*A[1][2] - A[0][2]*A[1][1];
35-
CeedScalar detA = A[0][0]*B11 + A[1][0]*B12 + A[2][0]*B13;
36-
37-
return detA;
38-
};
39-
#endif
40-
24+
#include "utils.h"
4125
// -----------------------------------------------------------------------------
4226
// This QFunction applies the mass operator for a vector field of 2 components.
4327
//
@@ -71,24 +55,17 @@ CEED_QFUNCTION(SetupMass3D)(void *ctx, CeedInt Q, const CeedScalar *const *in,
7155
const CeedScalar J[3][3] = {{dxdX[0][0][i], dxdX[1][0][i], dxdX[2][0][i]},
7256
{dxdX[0][1][i], dxdX[1][1][i], dxdX[2][1][i]},
7357
{dxdX[0][2][i], dxdX[1][2][i], dxdX[2][2][i]}};
74-
const CeedScalar detJ = ComputeDetMat(J);
75-
const CeedScalar u1[3] = {u[0][i], u[1][i], u[2][i]};
58+
const CeedScalar det_J = MatDet3x3(J);
59+
CeedScalar u1[3] = {u[0][i], u[1][i], u[2][i]}, v1[3];
7660
// *INDENT-ON*
7761
// Piola map: J^T*J*u*w/detJ
7862
// 1) Compute J^T * J
79-
CeedScalar JTJ[3][3];
80-
for (CeedInt j = 0; j < 3; j++) {
81-
for (CeedInt k = 0; k < 3; k++) {
82-
JTJ[j][k] = 0;
83-
for (CeedInt m = 0; m < 3; m++)
84-
JTJ[j][k] += J[m][j] * J[m][k];
85-
}
86-
}
63+
CeedScalar JT_J[3][3];
64+
AlphaMatTransposeMatMult3x3(1., J, J, JT_J);
8765
// 2) Compute J^T*J*u * w /detJ
66+
AlphaMatVecMult3x3(w[i]/det_J, JT_J, u1, v1);
8867
for (CeedInt k = 0; k < 3; k++) {
89-
v[k][i] = 0;
90-
for (CeedInt m = 0; m < 3; m++)
91-
v[k][i] += JTJ[k][m] * u1[m] * w[i]/detJ;
68+
v[k][i] = v1[k];
9269
}
9370
} // End of Quadrature Point Loop
9471

‎examples/Hdiv-mass/qfunctions/poisson-rhs2d.h‎

Lines changed: 3 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
#define POISSON_RHS2D_H
2222

2323
#include <math.h>
24+
#include "utils.h"
2425

2526
#ifndef M_PI
2627
#define M_PI 3.14159265358979323846
@@ -63,11 +64,8 @@ CEED_QFUNCTION(SetupRhs2D)(void *ctx, const CeedInt Q,
6364
CeedScalar ue[2] = {-M_PI*cos(M_PI*x) *sin(M_PI*y), -M_PI*sin(M_PI*x) *cos(M_PI*y)};
6465
//CeedScalar ue[2] = {x-y, x+y};
6566
CeedScalar rhs1[2];
66-
for (CeedInt k = 0; k < 2; k++) {
67-
rhs1[k] = 0;
68-
for (CeedInt m = 0; m < 2; m++)
69-
rhs1[k] += J[m][k] * ue[m];
70-
}
67+
AlphaMatTransposeVecMult2x2(1, J, ue, rhs1);
68+
7169
// Component 1
7270
true_soln[i+0*Q] = ue[0];
7371
rhs[i+0*Q] = rhs1[0] * w[i];

‎examples/Hdiv-mass/qfunctions/poisson-rhs3d.h‎

Lines changed: 2 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@
2121
#define POISSON_RHS3D_H
2222

2323
#include <math.h>
24-
24+
#include "utils.h"
2525
#ifndef M_PI
2626
#define M_PI 3.14159265358979323846
2727
#endif
@@ -66,11 +66,7 @@ CEED_QFUNCTION(SetupRhs3D)(void *ctx, const CeedInt Q,
6666
};
6767
//CeedScalar ue[3] = {x,y,z};
6868
CeedScalar rhs1[3];
69-
for (CeedInt k = 0; k < 3; k++) {
70-
rhs1[k] = 0;
71-
for (CeedInt m = 0; m < 3; m++)
72-
rhs1[k] += J[m][k] * ue[m];
73-
}
69+
AlphaMatTransposeVecMult3x3(1, J, ue, rhs1);
7470
// Component 1
7571
true_soln[i+0*Q] = ue[0];
7672
rhs[i+0*Q] = rhs1[0] * w[i];

0 commit comments

Comments
 (0)