FreeFem / FreeFem-sources

FreeFEM source code
https://freefem.org/
Other
775 stars 191 forks source link

Two questions about FreeFem++ #75

Closed zhaog6 closed 5 years ago

zhaog6 commented 5 years ago

Hello,

I'd like to ask if FreeFem++ support h- and p-adaptation respectively, and "trunc" in 3D(when I test Poission equations' convergence rate on norm H1 by using Th = trunc(Th, 1, split=2); to nestedly refine mesh in 3D, I find the convergence rate may be wrong) ? Thanks.

The code to test "trunc" in 3D is as following

load "msh3"
load "tetgen"
load "MUMPS"

macro fe P1 //

//-----------------+
//   Parameters    |
//-----------------+
int M = 32;
int nblayers = 1;

//-----------------------+
//    Exact Solution 1   |
//-----------------------+
func u = cos(2*pi*x) * cos(2*pi*y) * cos(2*pi*z);
func ux = - 2 * pi * sin(2*pi*x) * cos(2*pi*y) * cos(2*pi*z);
func uy = - 2 * pi * cos(2*pi*x) * sin(2*pi*y) * cos(2*pi*z);
func uz = - 2 * pi * cos(2*pi*x) * cos(2*pi*y) * sin(2*pi*z);
func uxx = - 4 * pi^2 * cos(2*pi*x) * cos(2*pi*y) * cos(2*pi*z);
func uyy = uxx;
func uzz = uxx;
func f = - uxx - uyy - uzz;
func g = u;

//-----------------+
//   Build Mesh    |
//-----------------+
border bd1(t = 0, 1) {x = t;   y = 0;   label=1;}
border bd2(t = 0, 1) {x = 1;   y = t;   label=2;}
border bd3(t = 0, 1) {x = 1-t; y = 1;   label=3;}
border bd4(t = 0, 1) {x = 0;   y = 1-t; label=4;}
mesh Th0 = buildmesh(bd1(M) + bd2(M) + bd3(M) + bd4(M));
mesh3 Th1 = movemesh23(Th0, transfo=[x, y, 0]);
mesh3 Th2 = movemesh23(Th0, transfo=[x, y, 1]);
mesh3 Th3 = movemesh23(Th0, transfo=[0, x, y]);
mesh3 Th4 = movemesh23(Th0, transfo=[1, x, y]);
mesh3 Th5 = movemesh23(Th0, transfo=[x, 0, y]);
mesh3 Th6 = movemesh23(Th0, transfo=[x, 1, y]);
mesh3[int] Th(nblayers);
Th[0] = Th1 + Th2 + Th3 + Th4 + Th5 + Th6;
Th[0] = tetg(Th[0]);
//plot(Th[0], cmm="Initial Mesh", wait=1);
for (int l = 1; l < nblayers; ++l)
{
    Th[l] = trunc(Th[l-1], 1, split=2);
}
cout << "Th[" << nblayers-1 << "].nt = " << Th[nblayers-1].nt << endl;
cout << "Th[" << nblayers-1 << "].nv = " << Th[nblayers-1].nv << endl;
//plot(Th[nblayers-1], cmm="Solve Mesh", wait=1);

//--------------------------------+
//   Build Finite Element Space   |
//--------------------------------+
fespace Vh(Th[nblayers-1], fe);
Vh uh = 0;

//-----------------------------+
//   Assemble Matrix and RHS   |
//-----------------------------+
varf vA(u, v) = int3d(Th[nblayers-1]) (dx(u)*dx(v) + dy(u)*dy(v) + dz(u)*dz(v)) + on(0, u=g);
varf vl(unused, v) = int3d(Th[nblayers-1]) (f * v) + on(0, unused=g);
matrix A = vA(Vh, Vh);
real[int] b = vl(0, Vh);

//-------------------+
//  Solve Equations  |
//-------------------+
set(A, solver=sparsesolver);
uh[] = A^-1 * b;

//-------------------------+
//   Compute Error Order   |
//-------------------------+
real H1Error = sqrt(int3d(Th[nblayers-1]) ((ux - dx(uh))^2 + (uy - dy(uh))^2 + (uz-dz(uh))^2 + (u - uh)^2));
cout << "=============================================" << endl;
cout << "H1Error = " << H1Error << endl;
cout << "=============================================" << endl;
//plot(uh, wait=1);

H1Error is 1.8384, but when changing code in line 11 into int nblayers = 2, H1Error is 1.12211. The convergence order would be bad (even though I refine grid again).

frederichecht commented 5 years ago

Hello, No, just h-adaptation.

Le 2 mars 2019 à 10:57, zhaog6 notifications@github.com a écrit :

Hello,

I'd like to ask if FreeFemm++ support h- and p-adaptation respectively, and "trunc" in 3d ? Thanks.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/FreeFem/FreeFem-sources/issues/75, or mute the thread https://github.com/notifications/unsubscribe-auth/AFb2jP3RX46QjuTXDKlYbuGVWr8Ihhb1ks5vSksfgaJpZM4baS_b.

zhaog6 commented 5 years ago

Thank you, Does the "trunc" function support nested refinement of 3D mesh?

prj- commented 5 years ago

Yes, through the optional parameter split.

zhaog6 commented 5 years ago

Thanks, I would like to ask another question. For h-adaptation, I find the mesh has moved when I use adaptmesh(...) to refine mesh according to "examples++-tutorial/AdaptResidualErrorIndicator.edp" (by given metric to refine). But I intend to refine mesh locally under fixed mesh to implement h-adaptation, which doesn't seem to have been mentioned in tutorial. So I would like to know which function to use.