Compare commits
2
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
 Jakub Cerveny
|
92e3a59b39 | ||
|
Jakub Cerveny
|
39efa26532 |
@@ -0,0 +1,108 @@
|
||||
#include "mfem.hpp"
|
||||
#include "error.hpp"
|
||||
|
||||
namespace mfem
|
||||
{
|
||||
|
||||
double CalculateH10Error2(GridFunction *sol, VectorCoefficient *exgrad,
|
||||
Array<double> *elemError, Array<int> *elemRef,
|
||||
int intOrder)
|
||||
{
|
||||
const FiniteElementSpace *fes = sol->FESpace();
|
||||
Mesh* mesh = fes->GetMesh();
|
||||
|
||||
Vector e_grad, a_grad, el_dofs, q_grad;
|
||||
DenseMatrix dshape, dshapet, Jinv;
|
||||
Array<int> vdofs;
|
||||
const FiniteElement *fe;
|
||||
ElementTransformation *transf;
|
||||
|
||||
int dim = mesh->Dimension();
|
||||
e_grad.SetSize(dim);
|
||||
a_grad.SetSize(dim);
|
||||
q_grad.SetSize(dim);
|
||||
Jinv.SetSize(dim);
|
||||
|
||||
double error = 0.0;
|
||||
if (elemError) { elemError->SetSize(mesh->GetNE()); }
|
||||
if (elemRef) { elemRef->SetSize(mesh->GetNE()); }
|
||||
|
||||
for (int i = 0; i < mesh->GetNE(); i++)
|
||||
{
|
||||
fe = fes->GetFE(i);
|
||||
int fdof = fe->GetDof();
|
||||
transf = mesh->GetElementTransformation(i);
|
||||
el_dofs.SetSize(fdof);
|
||||
dshape.SetSize(fdof, dim);
|
||||
dshapet.SetSize(fdof, dim);
|
||||
|
||||
fes->GetElementVDofs(i, vdofs);
|
||||
for (int k = 0; k < fdof; k++)
|
||||
{
|
||||
el_dofs(k) = (vdofs[k] >= 0) ? (*sol)(vdofs[k])
|
||||
: -(*sol)(-1-vdofs[k]);
|
||||
}
|
||||
|
||||
const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(), intOrder);
|
||||
|
||||
// integrate the H^1_0 error
|
||||
double el_err = 0.0, a_dxyz[3] = { 0, 0, 0 };
|
||||
for (int j = 0; j < ir.GetNPoints(); j++)
|
||||
{
|
||||
const IntegrationPoint &ip = ir.IntPoint(j);
|
||||
|
||||
transf->SetIntPoint(&ip);
|
||||
CalcInverse(transf->Jacobian(), Jinv);
|
||||
double w = ip.weight * transf->Weight();
|
||||
|
||||
exgrad->Eval(e_grad, *transf, ip);
|
||||
|
||||
fe->CalcDShape(ip, dshape);
|
||||
Mult(dshape, Jinv, dshapet);
|
||||
dshapet.MultTranspose(el_dofs, a_grad);
|
||||
|
||||
e_grad -= a_grad;
|
||||
el_err += w * (e_grad * e_grad);
|
||||
|
||||
// anisotropic indicators
|
||||
transf->Jacobian().MultTranspose(e_grad, q_grad);
|
||||
for (int k = 0; k < dim; k++)
|
||||
{
|
||||
a_dxyz[k] += w * (q_grad[k] * q_grad[k]);
|
||||
}
|
||||
}
|
||||
|
||||
error += el_err;
|
||||
if (elemError)
|
||||
{
|
||||
(*elemError)[i] = fabs(el_err);
|
||||
}
|
||||
|
||||
// determine what type of anisotropic refinement (if any) is suitable
|
||||
if (elemRef)
|
||||
{
|
||||
double sum = 0;
|
||||
for (int k = 0; k < dim; k++)
|
||||
{
|
||||
sum += a_dxyz[k];
|
||||
}
|
||||
|
||||
const double thresh = 0.2 * 3/dim;
|
||||
int ref = 0;
|
||||
for (int k = 0; k < dim; k++)
|
||||
{
|
||||
if (a_dxyz[k] / sum > thresh)
|
||||
{
|
||||
ref |= (1 << k);
|
||||
}
|
||||
}
|
||||
|
||||
(*elemRef)[i] = ref;
|
||||
}
|
||||
}
|
||||
|
||||
return error;
|
||||
}
|
||||
|
||||
|
||||
} // namespace mfem
|
||||
@@ -0,0 +1,26 @@
|
||||
#ifndef HPTEST_ERROR_HPP
|
||||
#define HPTEST_ERROR_HPP
|
||||
|
||||
namespace mfem
|
||||
{
|
||||
|
||||
/**
|
||||
* @brief Calcuate the square of the H^1_0 error of 'sol' against exact
|
||||
* solution 'exgrad'.
|
||||
*
|
||||
* @param sol Approximate solution.
|
||||
* @param exgrad Gradient of the exact solution.
|
||||
* @param elemError Optional array that receives per-element error.
|
||||
* @param elemRef Optional array that receives per-element anisotropy flag.
|
||||
* @param intOrder Integration rule order to use for evaluating the error.
|
||||
* @return
|
||||
*/
|
||||
double CalculateH10Error2(GridFunction *sol,
|
||||
VectorCoefficient *exgrad,
|
||||
Array<double> *elemError,
|
||||
Array<int> *elemRef,
|
||||
int intOrder);
|
||||
|
||||
} // namespace mfem
|
||||
|
||||
#endif // HPTEST_ERROR_HPP
|
||||
@@ -0,0 +1,250 @@
|
||||
#include "mfem.hpp"
|
||||
#include "exact.hpp"
|
||||
|
||||
namespace mfem
|
||||
{
|
||||
|
||||
// L-shape domain problem exact solution (2D)
|
||||
|
||||
double lshape_exsol(const Vector &p)
|
||||
{
|
||||
double x = p(0), y = p(1);
|
||||
double r = sqrt(x*x + y*y);
|
||||
double a = atan2(y, x);
|
||||
if (a < 0) a += 2*M_PI;
|
||||
return pow(r, 2.0/3.0) * sin(2.0*a/3.0);
|
||||
}
|
||||
|
||||
void lshape_exgrad(const Vector &p, Vector &grad)
|
||||
{
|
||||
double x = p(0), y = p(1);
|
||||
double a = atan2(y, x);
|
||||
if (a < 0) a += 2*M_PI;
|
||||
double theta23 = 2.0/3.0*a;
|
||||
double r23 = pow(x*x + y*y, 2.0/3.0);
|
||||
grad(0) = 2.0/3.0*x*sin(theta23)/(r23) - 2.0/3.0*y*cos(theta23)/(r23);
|
||||
grad(1) = 2.0/3.0*y*sin(theta23)/(r23) + 2.0/3.0*x*cos(theta23)/(r23);
|
||||
}
|
||||
|
||||
double lshape_laplace(const Vector &p)
|
||||
{
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
// Fichera's corner problem exact solution (3D)
|
||||
|
||||
double fichera_exsol(const Vector &p)
|
||||
{
|
||||
double x = p(0), y = p(1), z = p(2);
|
||||
return pow(x*x + y*y + z*z, 0.25);
|
||||
}
|
||||
|
||||
void fichera_exgrad(const Vector &p, Vector &grad)
|
||||
{
|
||||
double x = p(0), y = p(1), z = p(2);
|
||||
grad(0) = 0.5 * x * pow(x*x + y*y + z*z, -0.75);
|
||||
grad(1) = 0.5 * y * pow(x*x + y*y + z*z, -0.75);
|
||||
grad(2) = 0.5 * z * pow(x*x + y*y + z*z, -0.75);
|
||||
}
|
||||
|
||||
double fichera_laplace(const Vector &p)
|
||||
{
|
||||
double x = p(0), y = p(1), z = p(2);
|
||||
return -0.75 * pow(x*x + y*y + z*z, -0.75);
|
||||
}
|
||||
|
||||
|
||||
// inner layer problem exact solution (2D)
|
||||
|
||||
#if 1
|
||||
const double alpha = 200.0; // standard params
|
||||
const double center = -0.05;
|
||||
const double radius = 0.7;
|
||||
#elif 0
|
||||
const double alpha = 400.0; // cube centered
|
||||
const double center = 0.5;
|
||||
const double radius = 0.3;
|
||||
#elif 1
|
||||
const double alpha = 80.0; // nurbs ball
|
||||
const double center = -1;
|
||||
const double radius = 1.8;
|
||||
#else
|
||||
const double alpha = 80.0; // hcurl
|
||||
const double center = -0.05;
|
||||
const double radius = 0.7;
|
||||
#endif
|
||||
|
||||
template<typename T> T sqr(T x) { return x*x; }
|
||||
|
||||
double layer2_exsol(const Vector &p)
|
||||
{
|
||||
double x = p(0), y = p(1);
|
||||
double r = sqrt(sqr(x - center) + sqr(y - center));
|
||||
return atan(alpha * (r - radius));
|
||||
}
|
||||
|
||||
void layer2_exgrad(const Vector &p, Vector &grad)
|
||||
{
|
||||
double x = p(0), y = p(1);
|
||||
double r = sqrt(sqr(x - center) + sqr(y - center));
|
||||
double u = r * (sqr(alpha) * sqr(r - radius) + 1);
|
||||
grad(0) = alpha * (x - center) / u;
|
||||
grad(1) = alpha * (y - center) / u;
|
||||
}
|
||||
|
||||
double layer2_laplace(const Vector &p)
|
||||
{
|
||||
double x = p(0), y = p(1);
|
||||
double r = sqr(y - center) + sqr(x - center);
|
||||
double u = sqr(alpha) * sqr(sqrt(r) - radius) + 1;
|
||||
|
||||
return 2 * pow(alpha,3) * (sqrt(r) - radius) * sqr(y - center) / (r * sqr(u))
|
||||
+ alpha * sqr(y - center) / (pow(r, 1.5) * u)
|
||||
- 2 * alpha / (sqrt(r) * u)
|
||||
+ 2 * pow(alpha,3) * (sqrt(r) - radius) * sqr(x - center) / (r * sqr(u))
|
||||
+ alpha * sqr(x - center) / (pow(r, 1.5) * u);
|
||||
}
|
||||
|
||||
|
||||
// inner layer problem exact solution (3D)
|
||||
|
||||
double layer3_exsol(const Vector &p)
|
||||
{
|
||||
double x = p(0), y = p(1), z = p(2);
|
||||
double r = sqrt(sqr(x - center) + sqr(y - center) + sqr(z - center));
|
||||
return atan(alpha * (r - radius));
|
||||
}
|
||||
|
||||
void layer3_exgrad(const Vector &p, Vector &grad)
|
||||
{
|
||||
double x = p(0), y = p(1), z = p(2);
|
||||
double t1 = x * x;
|
||||
double t4 = center * center;
|
||||
double t6 = y * y;
|
||||
double t9 = z * z;
|
||||
double t13 = sqrt(t1 - 0.2e1 * x * center + 0.3e1 * t4 + t6
|
||||
- 0.2e1 * y * center + t9 - 0.2e1 * z * center);
|
||||
double t17 = alpha * alpha;
|
||||
double t19 = pow(t13 - radius, 0.2e1);
|
||||
grad(0) = alpha / t13 * (x - center) / (0.1e1 + t17 * t19);
|
||||
grad(1) = alpha / t13 * (y - center) / (0.1e1 + t17 * t19);
|
||||
grad(2) = alpha / t13 * (z - center) / (0.1e1 + t17 * t19);
|
||||
}
|
||||
|
||||
double layer3_laplace(const Vector &p)
|
||||
{
|
||||
double x = p(0), y = p(1), z = p(2);
|
||||
|
||||
double t1 = x * x;
|
||||
double t4 = center * center;
|
||||
double t6 = y * y;
|
||||
double t9 = z * z;
|
||||
double t12 = t1 - 0.2e1 * x * center + 0.3e1 * t4
|
||||
+ t6 - 0.2e1 * y * center + t9 - 0.2e1 * z * center;
|
||||
double t13 = sqrt(t12);
|
||||
double t16 = alpha / t13 / t12;
|
||||
double t18 = 0.4e1 * pow(x - center, 0.2e1);
|
||||
double t19 = alpha * alpha;
|
||||
double t20 = t13 - radius;
|
||||
double t21 = t20 * t20;
|
||||
double t23 = 0.1e1 + t19 * t21;
|
||||
double t24 = 0.1e1 / t23;
|
||||
double t34 = alpha * t19 / t12;
|
||||
double t35 = t23 * t23;
|
||||
double t36 = 0.1e1 / t35;
|
||||
double t42 = 0.4e1 * pow(y - center, 0.2e1);
|
||||
double t51 = 0.4e1 * pow(z - center, 0.2e1);
|
||||
double t59 = -t16 * t18 * t24 / 0.4e1
|
||||
+ 0.3e1 * alpha / t13 * t24
|
||||
- t34 * t18 * t36 * t20 / 0.2e1
|
||||
- t16 * t42 * t24 / 0.4e1
|
||||
- t34 * t42 * t36 * t20 / 0.2e1
|
||||
- t16 * t51 * t24 / 0.4e1
|
||||
- t34 * t51 * t36 * t20 / 0.2e1;
|
||||
|
||||
return -t59;
|
||||
}
|
||||
|
||||
|
||||
// hcurl
|
||||
|
||||
#if 0
|
||||
const double kappa = M_PI;
|
||||
|
||||
void hcurl_exsol(const Vector &x, Vector &E)
|
||||
{
|
||||
E(0) = sin(kappa * x(1));
|
||||
E(1) = sin(kappa * x(2));
|
||||
E(2) = sin(kappa * x(0));
|
||||
}
|
||||
|
||||
void hcurl_exrhs(const Vector &x, Vector &f)
|
||||
{
|
||||
f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
|
||||
f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
|
||||
f(2) = (1. + kappa * kappa) * sin(kappa * x(0));
|
||||
}
|
||||
#else
|
||||
|
||||
void hcurl_exsol(const Vector &x, Vector &E)
|
||||
{
|
||||
Vector a(x), b(x), c(x);
|
||||
a(0) = x(1); a(1) = x(2); // yz
|
||||
b(0) = x(0); b(1) = x(2); // xz
|
||||
|
||||
E(0) = layer2_exsol(a);
|
||||
E(1) = 0;//layer2_exsol(b);
|
||||
E(2) = 0;//layer2_exsol(c);
|
||||
}
|
||||
|
||||
void hcurl_exrhs(const Vector &x, Vector &f)
|
||||
{
|
||||
Vector a(x), b(x), c(x);
|
||||
a(0) = x(1); a(1) = x(2); // yz
|
||||
b(0) = x(0); b(1) = x(2); // xz
|
||||
|
||||
f(0) = layer2_laplace(a) + layer2_exsol(a);
|
||||
f(1) = 0;//layer2_laplace(b) + layer2_exsol(b);
|
||||
f(2) = 0;//layer2_laplace(c) + layer2_exsol(c);
|
||||
}
|
||||
|
||||
#endif
|
||||
|
||||
// hdiv
|
||||
|
||||
void hdiv_exsol(const Vector &p, Vector &F)
|
||||
{
|
||||
double x,y,z;
|
||||
int dim = p.Size();
|
||||
|
||||
x = p(0);
|
||||
y = p(1);
|
||||
if (dim == 3) { z = p(2); }
|
||||
|
||||
F(0) = cos(M_PI*x) * sin(M_PI*y);
|
||||
F(1) = cos(M_PI*y) * sin(M_PI*x);
|
||||
if (dim == 3) { F(2) = 0.0; }
|
||||
|
||||
(void) z;
|
||||
}
|
||||
|
||||
void hdiv_exrhs(const Vector &p, Vector &f)
|
||||
{
|
||||
double x,y,z;
|
||||
int dim = p.Size();
|
||||
|
||||
x = p(0);
|
||||
y = p(1);
|
||||
if (dim == 3) { z = p(2); }
|
||||
|
||||
double temp = 1 + 2*M_PI*M_PI;
|
||||
|
||||
f(0) = temp * cos(M_PI*x) * sin(M_PI*y);
|
||||
f(1) = temp * cos(M_PI*y) * sin(M_PI*x);
|
||||
if (dim == 3) { f(2) = 0; }
|
||||
|
||||
(void) z;
|
||||
}
|
||||
|
||||
} // namespace mfem
|
||||
@@ -0,0 +1,35 @@
|
||||
#ifndef HPTEST_EXACT_HPP
|
||||
#define HPTEST_EXACT_HPP
|
||||
|
||||
namespace mfem
|
||||
{
|
||||
|
||||
// classic L-shape domain problem with a singularity in reentrant corner
|
||||
double lshape_exsol(const Vector &p);
|
||||
void lshape_exgrad(const Vector &p, Vector &grad);
|
||||
double lshape_laplace(const Vector &p);
|
||||
|
||||
double fichera_exsol(const Vector &p);
|
||||
void fichera_exgrad(const Vector &p, Vector &grad);
|
||||
double fichera_laplace(const Vector &p);
|
||||
|
||||
// shock-like "inner layer" problem
|
||||
double layer2_exsol(const Vector &p);
|
||||
void layer2_exgrad(const Vector &p, Vector &grad);
|
||||
double layer2_laplace(const Vector &p);
|
||||
|
||||
// shock-like "inner layer" problem generalized to 3D
|
||||
double layer3_exsol(const Vector &p);
|
||||
void layer3_exgrad(const Vector &p, Vector &grad);
|
||||
double layer3_laplace(const Vector &p);
|
||||
|
||||
void hcurl_exsol(const Vector &x, Vector &E);
|
||||
void hcurl_exrhs(const Vector &x, Vector &f);
|
||||
|
||||
void hdiv_exsol(const Vector &x, Vector &E);
|
||||
void hdiv_exrhs(const Vector &x, Vector &f);
|
||||
|
||||
|
||||
} // namespace mfem
|
||||
|
||||
#endif // HPTEST_EXACT_HPP
|
||||
@@ -0,0 +1,82 @@
|
||||
MFEM mesh v1.0
|
||||
|
||||
#
|
||||
# MFEM Geometry Types (see mesh/geom.hpp):
|
||||
#
|
||||
# POINT = 0
|
||||
# SEGMENT = 1
|
||||
# TRIANGLE = 2
|
||||
# SQUARE = 3
|
||||
# TETRAHEDRON = 4
|
||||
# CUBE = 5
|
||||
#
|
||||
|
||||
dimension
|
||||
3
|
||||
|
||||
elements
|
||||
7
|
||||
1 5 0 1 4 3 9 10 13 12
|
||||
1 5 1 2 5 4 10 11 14 13
|
||||
1 5 3 4 7 6 12 13 16 15
|
||||
1 5 4 5 8 7 13 14 17 16
|
||||
1 5 9 10 13 12 18 19 21 20
|
||||
1 5 12 13 16 15 20 21 24 23
|
||||
1 5 13 14 17 16 21 22 25 24
|
||||
|
||||
boundary
|
||||
24
|
||||
1 3 0 3 12 9
|
||||
1 3 3 6 15 12
|
||||
1 3 9 12 20 18
|
||||
1 3 12 15 23 20
|
||||
1 3 2 5 14 11
|
||||
1 3 5 8 17 14
|
||||
1 3 10 13 21 19
|
||||
1 3 14 17 25 22
|
||||
1 3 0 1 10 9
|
||||
1 3 1 2 11 10
|
||||
1 3 9 10 19 18
|
||||
1 3 13 14 22 21
|
||||
1 3 6 7 16 15
|
||||
1 3 7 8 17 16
|
||||
1 3 15 16 24 23
|
||||
1 3 16 17 25 24
|
||||
1 3 0 1 4 3
|
||||
1 3 1 2 5 4
|
||||
1 3 3 4 7 6
|
||||
1 3 4 5 8 7
|
||||
1 3 18 19 21 20
|
||||
1 3 10 11 14 13
|
||||
1 3 20 21 24 23
|
||||
1 3 21 22 25 24
|
||||
|
||||
vertices
|
||||
26
|
||||
3
|
||||
-1 -1 -1
|
||||
0 -1 -1
|
||||
1 -1 -1
|
||||
-1 0 -1
|
||||
0 0 -1
|
||||
1 0 -1
|
||||
-1 1 -1
|
||||
0 1 -1
|
||||
1 1 -1
|
||||
-1 -1 0
|
||||
0 -1 0
|
||||
1 -1 0
|
||||
-1 0 0
|
||||
0 0 0
|
||||
1 0 0
|
||||
-1 1 0
|
||||
0 1 0
|
||||
1 1 0
|
||||
-1 -1 1
|
||||
0 -1 1
|
||||
-1 0 1
|
||||
0 0 1
|
||||
1 0 1
|
||||
-1 1 1
|
||||
0 1 1
|
||||
1 1 1
|
||||
@@ -0,0 +1,614 @@
|
||||
// hp-Refinement Demo
|
||||
//
|
||||
// Compile with: make hptest
|
||||
//
|
||||
//
|
||||
// Description: This is a demo of the hp-refinement capability of MFEM.
|
||||
// Two benchmark problems with a known exact solution are
|
||||
// solved on a sequence of meshes where both the size (h) and the
|
||||
// polynomial order (p) of elements is adapted.
|
||||
//
|
||||
|
||||
#include "mfem.hpp"
|
||||
#include <fstream>
|
||||
|
||||
#include "exact.hpp"
|
||||
#include "util.hpp"
|
||||
#include "error.hpp"
|
||||
|
||||
using namespace std;
|
||||
using namespace mfem;
|
||||
|
||||
const char* keys = "Rjlmc*******";
|
||||
|
||||
struct HPRefinement : public Refinement
|
||||
{
|
||||
int orders[4];
|
||||
|
||||
HPRefinement() = default;
|
||||
|
||||
HPRefinement(int index, int type = 7)
|
||||
: Refinement(index, type)
|
||||
{
|
||||
orders[0] = orders[1] = orders[2] = orders[3] = 0;
|
||||
|
||||
}
|
||||
};
|
||||
|
||||
struct HPCandidate
|
||||
{
|
||||
double err;
|
||||
int dof;
|
||||
int orders[4];
|
||||
};
|
||||
|
||||
|
||||
void MakeConforming(GridFunction &sol)
|
||||
{
|
||||
FiniteElementSpace* fes = sol.FESpace();
|
||||
const SparseMatrix* P = fes->GetConformingProlongation();
|
||||
const SparseMatrix* Q = fes->GetConformingRestrictionInterpolation();
|
||||
if (P)
|
||||
{
|
||||
Vector X;
|
||||
X.SetSize(Q->Height());
|
||||
Q->Mult(sol, X);
|
||||
P->Mult(X, sol);
|
||||
}
|
||||
}
|
||||
|
||||
int HalfOrder(int p)
|
||||
{
|
||||
return max(1, p/2);
|
||||
}
|
||||
|
||||
bool ContainsVertex(Mesh *mesh, int elem, const Vertex& vert)
|
||||
{
|
||||
Array<int> v;
|
||||
mesh->GetElementVertices(elem, v);
|
||||
for (int j = 0; j < v.Size(); j++)
|
||||
{
|
||||
double* vertex = mesh->GetVertex(v[j]);
|
||||
double dist = 0.0;
|
||||
for (int l = 0; l < 2; l++)
|
||||
{
|
||||
double d = vert(l) - vertex[l];
|
||||
dist += d*d;
|
||||
}
|
||||
if (dist == 0) { return true; }
|
||||
}
|
||||
return false;
|
||||
}
|
||||
|
||||
void Solve(FiniteElementSpace *fespace, GridFunction *sln, Coefficient *exsol, Coefficient *rhs,
|
||||
bool pa, int int_order, bool relaxed_hp)
|
||||
{
|
||||
DomainLFIntegrator *dlfi = new DomainLFIntegrator(*rhs);
|
||||
Geometry::Type geom = fespace->GetMesh()->GetElementGeometry(0);
|
||||
dlfi->SetIntRule(&IntRules.Get(geom, int_order));
|
||||
|
||||
// Assemble the linear form. The right hand side is manufactured
|
||||
// so that the solution is the analytic solution.
|
||||
LinearForm lf(fespace);
|
||||
lf.AddDomainIntegrator(dlfi);
|
||||
lf.Assemble();
|
||||
|
||||
double sigma = -1;
|
||||
double kappa = 1;
|
||||
|
||||
// Assemble the bilinear form.
|
||||
BilinearForm bf(fespace);
|
||||
if (pa) { bf.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
|
||||
bf.AddDomainIntegrator(new DiffusionIntegrator());
|
||||
if (relaxed_hp) {
|
||||
bf.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(sigma, kappa));
|
||||
}
|
||||
bf.Assemble();
|
||||
|
||||
// Set Dirichlet boundary values in the GridFunction x.
|
||||
// Determine the list of Dirichlet true DOFs in the linear system.
|
||||
Array<int> ess_bdr(fespace->GetMesh()->bdr_attributes.Max());
|
||||
ess_bdr = 1;
|
||||
*sln = 0; // FIXME
|
||||
sln->ProjectBdrCoefficient(*exsol, ess_bdr);
|
||||
Array<int> ess_tdof_list;
|
||||
fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
|
||||
|
||||
// 16. Create the linear system: eliminate boundary conditions, constrain
|
||||
// hanging nodes and possibly apply other transformations. The system
|
||||
// will be solved for true (unconstrained) DOFs only.
|
||||
OperatorPtr A;
|
||||
Vector B, X;
|
||||
|
||||
const int copy_interior = 1;
|
||||
bf.FormLinearSystem(ess_tdof_list, *sln, lf, A, X, B, copy_interior);
|
||||
|
||||
// 17. Solve the linear system A X = B.
|
||||
if (!pa)
|
||||
{
|
||||
#ifndef MFEM_USE_SUITESPARSE
|
||||
// Use a simple symmetric Gauss-Seidel preconditioner with PCG.
|
||||
GSSmoother M((SparseMatrix&)(*A));
|
||||
PCG(*A, M, B, X, 3, 2000, 1e-30, 0.0);
|
||||
#else
|
||||
// If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
|
||||
UMFPackSolver umf_solver;
|
||||
umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
|
||||
umf_solver.SetOperator(*A);
|
||||
umf_solver.Mult(B, X);
|
||||
#endif
|
||||
}
|
||||
else // No preconditioning for now in partial assembly mode.
|
||||
{
|
||||
CG(*A, B, X, 3, 2000, 1e-12, 0.0);
|
||||
}
|
||||
|
||||
// 18. After solving the linear system, reconstruct the solution as a
|
||||
// finite element GridFunction. Constrained nodes are interpolated
|
||||
// from true DOFs (it may therefore happen that x.Size() >= X.Size()).
|
||||
bf.RecoverFEMSolution(X, lf, *sln);
|
||||
}
|
||||
|
||||
|
||||
struct Solution
|
||||
{
|
||||
Mesh mesh;
|
||||
FiniteElementSpace fes;
|
||||
GridFunction sol;
|
||||
Array<double> elemError;
|
||||
|
||||
Solution(FiniteElementSpace &fespace, bool h_refined, int order_increase)
|
||||
: mesh(*(fespace.GetMesh())), fes(fespace, &mesh), sol(&fes)
|
||||
{
|
||||
Array<Refinement> mesh_refinements;
|
||||
for (int i = 0; i < fespace.GetNE(); i++)
|
||||
{
|
||||
int o, p = fespace.GetElementOrder(i);
|
||||
if (h_refined)
|
||||
{
|
||||
mesh_refinements.Append(Refinement(i));
|
||||
o = min(p, HalfOrder(p) + order_increase);
|
||||
}
|
||||
else
|
||||
{
|
||||
o = p + order_increase;
|
||||
}
|
||||
fes.SetElementOrder(i, o);
|
||||
}
|
||||
fes.Update(false);
|
||||
sol.Update();
|
||||
|
||||
if (h_refined)
|
||||
{
|
||||
mesh.GeneralRefinement(mesh_refinements, -1, 2);
|
||||
fes.Update(false);
|
||||
sol.Update();
|
||||
}
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
void FindHPRef(int elem, FiniteElementSpace *fes, Array<double> elemError, int n,
|
||||
Array<Solution*> solution, int max_order, std::map<int, HPRefinement> *hp_refs)
|
||||
{
|
||||
Mesh* mesh_h = solution[1]->fes.GetMesh();
|
||||
|
||||
// Find sons of elem
|
||||
int sons[4];
|
||||
const CoarseFineTransformations tr = mesh_h->GetRefinementTransforms();
|
||||
int l = 0;
|
||||
for (int i = 0; i < mesh_h->GetNE(); i++)
|
||||
{
|
||||
int j = tr.embeddings[i].parent;
|
||||
if (j == elem)
|
||||
{
|
||||
sons[l] = i;
|
||||
l++;
|
||||
}
|
||||
}
|
||||
|
||||
int o = fes->GetElementOrder(elem);
|
||||
int op = solution[0]->fes.GetElementOrder(elem);
|
||||
|
||||
double s_err = sqrt(elemError[elem]);
|
||||
cout << "Element " << elem << " (order " << o << "): err = " << s_err << ", dof = " << o*o << "\n ";
|
||||
|
||||
|
||||
// initialize candidates
|
||||
int n_cand = (n-1)*(n-1)*(n-1)*(n-1) + 1;
|
||||
HPCandidate candidate[n_cand];
|
||||
for (int id = 0; id < n_cand; id++)
|
||||
{
|
||||
candidate[id].err = 0;
|
||||
candidate[id].dof = 0;
|
||||
}
|
||||
|
||||
// p-candidate:
|
||||
candidate[0].err = solution[0]->elemError[elem];
|
||||
candidate[0].dof = op*op;
|
||||
candidate[0].orders[0] = op;
|
||||
|
||||
// hp-candidates
|
||||
int cand_id = 1;
|
||||
int k[4];
|
||||
for (k[0] = 1; k[0] < n; k[0]++)
|
||||
{
|
||||
for (k[1] = 1; k[1] < n; k[1]++)
|
||||
{
|
||||
for (k[2] = 1; k[2] < n; k[2]++)
|
||||
{
|
||||
for (k[3] = 1; k[3] < n; k[3]++)
|
||||
{
|
||||
for (int son = 0; son < 4; son++)
|
||||
{
|
||||
int oh = solution[k[son]]->fes.GetElementOrder(sons[son]);
|
||||
candidate[cand_id].err += solution[k[son]]->elemError[sons[son]];
|
||||
candidate[cand_id].dof += oh*oh;
|
||||
candidate[cand_id].orders[son] = oh;
|
||||
}
|
||||
cand_id++;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double max_rate = -1000.0;
|
||||
int best_id = 1;
|
||||
|
||||
for (int id = 0; id < n_cand; id++)
|
||||
{
|
||||
// define rate between error decrease and DOFs increase
|
||||
double rate = (s_err - sqrt(candidate[id].err)) / (candidate[id].dof - o*o);
|
||||
// double rate = (elemError[elem] - (candidate[id].err)) / (candidate[id].dof - o*o);
|
||||
|
||||
// // print all candidates
|
||||
// if (id == 0)
|
||||
// {
|
||||
// cout << "Candidate: " << id << ", err = " << sqrt(candidate[id].err) << ", rate = " << rate << ", dof = " << candidate[id].dof
|
||||
// << ", orders = " << candidate[id].orders[0] << "\n ";
|
||||
// }
|
||||
// else
|
||||
// {
|
||||
// cout << "Candidate: " << id << ", err = " << sqrt(candidate[id].err) << ", rate = " << rate << ", dof = " << candidate[id].dof
|
||||
// << ", orders = " << candidate[id].orders[0] << " " << candidate[id].orders[1]
|
||||
// << " " << candidate[id].orders[2] << " " << candidate[id].orders[3] << "\n ";
|
||||
// }
|
||||
|
||||
// throw away candidates with no error decrease or no DOF increase
|
||||
if ((elemError[elem] < (candidate[id].err)) || (candidate[id].dof <= o*o))
|
||||
continue;
|
||||
|
||||
// find candidate with highest rate
|
||||
if (rate > max_rate && candidate[id].orders[0] <= max_order)
|
||||
{
|
||||
max_rate = rate;
|
||||
best_id = id;
|
||||
}
|
||||
}
|
||||
|
||||
// print the best candidate
|
||||
if (best_id == 0)
|
||||
{
|
||||
cout << "Best candidate: " << best_id << ", err = " << sqrt(candidate[best_id].err) << ", dof = "
|
||||
<< candidate[best_id].dof << ", orders = " << candidate[best_id].orders[0] << "\n ";
|
||||
}
|
||||
else
|
||||
{
|
||||
cout << "Best candidate: " << best_id << ", err = " << sqrt(candidate[best_id].err) << ", dof = "
|
||||
<< candidate[best_id].dof << ", orders = " << candidate[best_id].orders[0] << " " << candidate[best_id].orders[1]
|
||||
<< " " << candidate[best_id].orders[2] << " " << candidate[best_id].orders[3] << "\n ";
|
||||
}
|
||||
|
||||
// Put the best candidate into hp_refinements
|
||||
(*hp_refs)[elem].index = elem;
|
||||
if (best_id > 0) { (*hp_refs)[elem].ref_type = 7; }
|
||||
else { (*hp_refs)[elem].ref_type = 0; }
|
||||
for (int i = 0; i < 4; i++)
|
||||
{
|
||||
(*hp_refs)[elem].orders[i] = candidate[best_id].orders[i];
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
// Parse command-line options.
|
||||
int dim = 2;
|
||||
int problem = 1;
|
||||
int order = 1;
|
||||
double ref_threshold = 0.7;
|
||||
bool aniso = false;
|
||||
bool hp = true;
|
||||
int n_enriched = 4;
|
||||
int max_order = 12;
|
||||
int int_order = 10;
|
||||
bool relaxed_hp = false;
|
||||
const char *conv_file = "conv.err";
|
||||
bool wait = false;
|
||||
int nc_limit = 2;
|
||||
bool pa = false;
|
||||
const char *device_config = "cpu";
|
||||
bool visualization = true;
|
||||
|
||||
OptionsParser args(argc, argv);
|
||||
args.AddOption(&problem, "-p", "--problem",
|
||||
"Problem type: 0 = L-shaped, 1 = inner layer.");
|
||||
args.AddOption(&dim, "-dim", "--dimension", "Dimension (2 or 3).");
|
||||
args.AddOption(&order, "-o", "--order",
|
||||
"Initial mesh finite element order (polynomial degree).");
|
||||
args.AddOption(&hp, "-hp", "--hp", "-no-hp", "--no-hp",
|
||||
"Enable hp refinement.");
|
||||
args.AddOption(&relaxed_hp, "-x", "--relaxed-hp", "-no-x", "--no-relaxed-hp",
|
||||
"Set relaxed hp conformity.");
|
||||
args.AddOption(&n_enriched, "-n", "--n_enriched",
|
||||
"Set number of enriched spaces (minimal value = 2).");
|
||||
args.AddOption(&conv_file, "-f", "--file",
|
||||
"Convergence file to use.");
|
||||
args.AddOption(&ref_threshold, "-rt", "--ref-threshold",
|
||||
"Refine elements with error larger than threshold * max_error.");
|
||||
args.AddOption(&nc_limit, "-nc", "--nc-limit",
|
||||
"Set maximum difference of refinement levels of adjacent elements.");
|
||||
args.AddOption(&wait, "-w", "--wait", "-no-w", "--no-wait",
|
||||
"Wait for user input after each iteration.");
|
||||
args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
|
||||
"--no-partial-assembly", "Enable Partial Assembly.");
|
||||
args.AddOption(&device_config, "-dev", "--device",
|
||||
"Device configuration string, see Device::Configure().");
|
||||
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
|
||||
"--no-visualization",
|
||||
"Enable or disable GLVis visualization.");
|
||||
args.Parse();
|
||||
if (!args.Good())
|
||||
{
|
||||
args.PrintUsage(cout);
|
||||
return 1;
|
||||
}
|
||||
args.PrintOptions(cout);
|
||||
|
||||
MFEM_VERIFY(dim >= 2 && dim <= 3, "Invalid dimension.");
|
||||
MFEM_VERIFY(problem >= 0 && problem <= 1, "Invalid problem type.");
|
||||
MFEM_VERIFY(n_enriched >= 2, "Invalid number of enriched spaces.");
|
||||
|
||||
// Enable hardware devices such as GPUs, and programming models such as
|
||||
// CUDA, OCCA, RAJA and OpenMP based on command line options.
|
||||
Device device(device_config);
|
||||
device.Print();
|
||||
|
||||
// Load and adjust the Mesh
|
||||
const char *mesh_file =
|
||||
problem ? ((dim == 3) ? "layer-hex.mesh" : "layer-quad.mesh")
|
||||
: ((dim == 3) ? "fichera-hex.mesh" : "lshape-quad.mesh");
|
||||
|
||||
Mesh mesh(mesh_file, 1, 1);
|
||||
|
||||
if (mesh.NURBSext)
|
||||
{
|
||||
for (int i = 0; i < 2; i++)
|
||||
{
|
||||
mesh.UniformRefinement();
|
||||
}
|
||||
mesh.SetCurvature(2);
|
||||
}
|
||||
mesh.EnsureNCMesh(true);
|
||||
|
||||
mesh.UniformRefinement();
|
||||
|
||||
// We don't support mixed meshes at the moment
|
||||
MFEM_VERIFY(mesh.GetNumGeometries(dim) == 1, "Mixed meshes not supported.");
|
||||
//Geometry::Type geom = mesh.GetElementGeometry(0);
|
||||
|
||||
// Prepare exact solution Coefficients
|
||||
FunctionCoefficient exsol(
|
||||
problem ? ((dim == 3) ? layer3_exsol : layer2_exsol)
|
||||
: ((dim == 3) ? fichera_exsol : lshape_exsol) );
|
||||
|
||||
VectorFunctionCoefficient exgrad(dim,
|
||||
problem ? ((dim == 3) ? layer3_exgrad : layer2_exgrad)
|
||||
: ((dim == 3) ? fichera_exgrad : lshape_exgrad) );
|
||||
|
||||
FunctionCoefficient rhs(
|
||||
problem ? ((dim == 3) ? layer3_laplace : layer2_laplace)
|
||||
: ((dim == 3) ? fichera_laplace : lshape_laplace) );
|
||||
|
||||
// Define a finite element space on the mesh. Initially the polynomial
|
||||
// order is constant everywhere.
|
||||
H1_FECollection fec(order, dim);
|
||||
FiniteElementSpace fespace(&mesh, &fec);
|
||||
fespace.SetRelaxedHpConformity(relaxed_hp);
|
||||
|
||||
// All boundary attributes will be used for essential (Dirichlet) BC.
|
||||
MFEM_VERIFY(mesh.bdr_attributes.Size() > 0,
|
||||
"Boundary attributes required in the mesh.");
|
||||
|
||||
// Connect to GLVis.
|
||||
socketstream sol_sock, ord_sock, dbg_sock[3], err_sock;
|
||||
|
||||
std::ofstream conv(conv_file);
|
||||
|
||||
// The main AMR loop. In each iteration we solve the problem on the
|
||||
// current mesh, visualize the solution, and refine the mesh.
|
||||
const int max_dofs = 100000;
|
||||
for (int it = 0; ; it++)
|
||||
{
|
||||
int cdofs = fespace.GetTrueVSize();
|
||||
cout << "\nAMR iteration " << it << endl;
|
||||
cout << "Number of unknowns: " << cdofs << endl;
|
||||
|
||||
// Solve for the current mesh: (h, p)
|
||||
GridFunction sol(&fespace);
|
||||
Solve(&fespace, &sol, &exsol, &rhs, pa, int_order, relaxed_hp);
|
||||
|
||||
// Calculate the H^1_0 errors of elements as well as the total error.
|
||||
Array<int> ref_type;
|
||||
Array<double> elemError;
|
||||
double error = sqrt(CalculateH10Error2(&sol, &exgrad, &elemError, &ref_type, int_order));
|
||||
double err_max = sqrt(elemError.Max());
|
||||
|
||||
// 19. Send solution by socket to the GLVis server.
|
||||
if (visualization)
|
||||
{
|
||||
GridFunction *vis_x = ProlongToMaxOrder(&sol);
|
||||
VisualizeField(sol_sock, *vis_x, "Solution", keys, 600, 500, 0, 70);
|
||||
delete vis_x;
|
||||
|
||||
GridFunction projsol(&fespace);
|
||||
Vector tmp = sol;
|
||||
projsol.ProjectCoefficient(exsol);
|
||||
MakeConforming(projsol);
|
||||
projsol -= tmp;
|
||||
vis_x = ProlongToMaxOrder(&projsol);
|
||||
VisualizeField(err_sock, *vis_x, "Error (Projection)", keys, 600, 500, 0, 70);
|
||||
delete vis_x;
|
||||
|
||||
L2_FECollection l2fec(0, dim);
|
||||
FiniteElementSpace l2fes(&mesh, &l2fec);
|
||||
GridFunction orders(&l2fes);
|
||||
for (int i = 0; i < orders.Size(); i++)
|
||||
{
|
||||
orders(i) = fespace.GetElementOrder(i);
|
||||
}
|
||||
VisualizeField(ord_sock, orders, "Orders", keys, 600, 500, 0, 620);
|
||||
//ord_sock << "valuerange 1 5\n" << flush;
|
||||
}
|
||||
|
||||
if (cdofs > max_dofs)
|
||||
{
|
||||
cout << "Reached the maximum number of dofs. Stop." << endl;
|
||||
break;
|
||||
}
|
||||
|
||||
// Save dofs and error for convergence plot
|
||||
conv << cdofs << " " << error << endl;
|
||||
cout << cdofs << " " << error << endl;
|
||||
|
||||
// Project the exact solution to h-refined and p-refined versions of the
|
||||
// mesh and determine whether to refine elements in 'h' or in 'p'.
|
||||
if (hp)
|
||||
{
|
||||
// Prepare enriched spaces
|
||||
Array<Solution*> solution(n_enriched);
|
||||
solution[0] = new Solution(fespace, false, 1); // p-refined space
|
||||
for (int k = 0; k < n_enriched - 1; k++)
|
||||
{
|
||||
solution[k+1] = new Solution(fespace, true, k); // h-refined spaces
|
||||
}
|
||||
// Solve for solution on enriched spaces
|
||||
for (int k = 0; k < n_enriched; k++)
|
||||
{
|
||||
Solve(&(solution[k]->fes), &(solution[k]->sol), &exsol, &rhs, pa, int_order, relaxed_hp);
|
||||
CalculateH10Error2(&(solution[k]->sol), &exgrad, &(solution[k]->elemError), &ref_type, int_order);
|
||||
}
|
||||
|
||||
int h_refined = 0, p_refined = 0;
|
||||
|
||||
Array<Refinement> refinements;
|
||||
std::map<int, HPRefinement> hp_refs;
|
||||
|
||||
for (int i = 0; i < mesh.GetNE(); i++)
|
||||
{
|
||||
if (sqrt(elemError[i]) > ref_threshold * err_max)
|
||||
{
|
||||
// Find the best hp refinement
|
||||
FindHPRef(i, &fespace, elemError, n_enriched, solution, max_order, &hp_refs);
|
||||
|
||||
if (hp_refs[i].ref_type > 0)
|
||||
{
|
||||
refinements.Append(Refinement(i));
|
||||
h_refined++;
|
||||
cout << "=> h-refined" << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
int p = fespace.GetElementOrder(i);
|
||||
fespace.SetElementOrder(i, p+1);
|
||||
p_refined++;
|
||||
cout << "=> p-refined" << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Update the space, interpolate the solution. FIXME
|
||||
fespace.Update(false);
|
||||
|
||||
// h-refine elements
|
||||
mesh.GeneralRefinement(refinements, -1, nc_limit);
|
||||
fespace.Update(false);
|
||||
|
||||
if (refinements.Size())
|
||||
{
|
||||
// Assign sons to all parents
|
||||
std::map<int, std::vector<int>> ref_sons;
|
||||
const CoarseFineTransformations tr = mesh.GetRefinementTransforms();
|
||||
for (int i = 0; i < mesh.GetNE(); i++)
|
||||
{
|
||||
int j = tr.embeddings[i].parent;
|
||||
ref_sons[j].push_back(i);
|
||||
}
|
||||
|
||||
// set orders for h-refined elements
|
||||
for (int i = 0; i < mesh.GetNE(); i++)
|
||||
{
|
||||
int j = tr.embeddings[i].parent;
|
||||
if (sqrt(elemError[j]) > ref_threshold * err_max)
|
||||
{
|
||||
if (hp_refs[j].ref_type > 0)
|
||||
{
|
||||
for (int k = 0; k < 4; k++)
|
||||
{
|
||||
fespace.SetElementOrder(ref_sons[j][k], hp_refs[j].orders[k]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
cout << "\nh-refined = " << h_refined
|
||||
<< ", p-refined = " << p_refined << endl;
|
||||
|
||||
for (int k = 0; k < n_enriched; k++)
|
||||
{
|
||||
delete solution[k];
|
||||
}
|
||||
|
||||
if (wait)
|
||||
{
|
||||
cout << "Press ENTER to continue...";
|
||||
cin.get();
|
||||
}
|
||||
|
||||
}
|
||||
else // !hp
|
||||
{
|
||||
if (wait)
|
||||
{
|
||||
cout << "Press ENTER to continue...";
|
||||
cin.get();
|
||||
}
|
||||
|
||||
Array<Refinement> refinements;
|
||||
double err_max = sqrt(elemError.Max());
|
||||
for (int i = 0; i < mesh.GetNE(); i++)
|
||||
{
|
||||
if (sqrt(elemError[i]) > ref_threshold * err_max)
|
||||
{
|
||||
int type = aniso ? ref_type[i] : 7;
|
||||
refinements.Append(Refinement(i, type));
|
||||
}
|
||||
}
|
||||
mesh.GeneralRefinement(refinements, -1, nc_limit);
|
||||
}
|
||||
|
||||
// Update the space, interpolate the solution.
|
||||
fespace.Update(false);
|
||||
|
||||
sol.Update();
|
||||
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,9 @@
|
||||
MFEM INLINE mesh v1.0
|
||||
|
||||
type = hex
|
||||
nx = 1
|
||||
ny = 1
|
||||
nz = 1
|
||||
sx = 1.0
|
||||
sy = 1.0
|
||||
sz = 1.0
|
||||
@@ -0,0 +1,34 @@
|
||||
MFEM mesh v1.0
|
||||
|
||||
#
|
||||
# MFEM Geometry Types (see mesh/geom.hpp):
|
||||
#
|
||||
# POINT = 0
|
||||
# SEGMENT = 1
|
||||
# TRIANGLE = 2
|
||||
# SQUARE = 3
|
||||
# TETRAHEDRON = 4
|
||||
# CUBE = 5
|
||||
#
|
||||
|
||||
dimension
|
||||
2
|
||||
|
||||
elements
|
||||
1
|
||||
1 3 0 1 2 3
|
||||
|
||||
boundary
|
||||
4
|
||||
1 1 0 1
|
||||
1 1 1 2
|
||||
1 1 2 3
|
||||
1 1 3 0
|
||||
|
||||
vertices
|
||||
4
|
||||
2
|
||||
0 0
|
||||
1 0
|
||||
1 1
|
||||
0 1
|
||||
@@ -0,0 +1,50 @@
|
||||
MFEM mesh v1.0
|
||||
|
||||
#
|
||||
# MFEM Geometry Types (see mesh/geom.hpp):
|
||||
#
|
||||
# POINT = 0
|
||||
# SEGMENT = 1
|
||||
# TRIANGLE = 2
|
||||
# SQUARE = 3
|
||||
# TETRAHEDRON = 4
|
||||
# CUBE = 5
|
||||
#
|
||||
|
||||
dimension
|
||||
2
|
||||
|
||||
elements
|
||||
8
|
||||
1 2 1 8 0
|
||||
1 2 1 2 8
|
||||
1 2 3 8 2
|
||||
1 2 3 4 8
|
||||
1 2 5 8 4
|
||||
1 2 5 6 8
|
||||
1 2 7 8 6
|
||||
1 2 7 0 8
|
||||
|
||||
boundary
|
||||
8
|
||||
1 1 0 1
|
||||
1 1 1 2
|
||||
1 1 2 3
|
||||
1 1 3 4
|
||||
1 1 4 5
|
||||
1 1 5 6
|
||||
1 1 6 7
|
||||
1 1 7 0
|
||||
|
||||
vertices
|
||||
9
|
||||
2
|
||||
0 0
|
||||
0.5 0
|
||||
1 0
|
||||
1 0.5
|
||||
1 1
|
||||
0.5 1
|
||||
0 1
|
||||
0 0.5
|
||||
0.5 0.5
|
||||
@@ -0,0 +1,44 @@
|
||||
MFEM mesh v1.0
|
||||
|
||||
#
|
||||
# MFEM Geometry Types (see mesh/geom.hpp):
|
||||
#
|
||||
# POINT = 0
|
||||
# SEGMENT = 1
|
||||
# TRIANGLE = 2
|
||||
# SQUARE = 3
|
||||
# TETRAHEDRON = 4
|
||||
# CUBE = 5
|
||||
#
|
||||
|
||||
dimension
|
||||
2
|
||||
|
||||
elements
|
||||
3
|
||||
1 3 0 1 2 3
|
||||
1 3 5 0 3 4
|
||||
1 3 6 7 0 5
|
||||
|
||||
boundary
|
||||
8
|
||||
1 1 0 1
|
||||
1 1 1 2
|
||||
1 1 2 3
|
||||
1 1 3 4
|
||||
1 1 4 5
|
||||
1 1 5 6
|
||||
1 1 6 7
|
||||
1 1 7 0
|
||||
|
||||
vertices
|
||||
8
|
||||
2
|
||||
0 0
|
||||
1 0
|
||||
1 1
|
||||
0 1
|
||||
-1 1
|
||||
-1 0
|
||||
-1 -1
|
||||
0 -1
|
||||
@@ -0,0 +1,66 @@
|
||||
# Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
|
||||
# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
|
||||
# LICENSE and NOTICE for details. LLNL-CODE-806117.
|
||||
#
|
||||
# This file is part of the MFEM library. For more information and source code
|
||||
# availability visit https://mfem.org.
|
||||
#
|
||||
# MFEM is free software; you can redistribute it and/or modify it under the
|
||||
# terms of the BSD-3 license. We welcome feedback and contributions, see file
|
||||
# CONTRIBUTING.md for details.
|
||||
|
||||
# Use the MFEM build directory
|
||||
MFEM_DIR ?= ../..
|
||||
MFEM_BUILD_DIR ?= ../..
|
||||
SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/hpfem/,)
|
||||
CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
|
||||
# Use the MFEM install directory
|
||||
# MFEM_INSTALL_DIR = ../../mfem
|
||||
# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
|
||||
|
||||
MFEM_LIB_FILE = mfem_is_not_built
|
||||
-include $(CONFIG_MK)
|
||||
|
||||
SEQ_MINIAPPS = hptest
|
||||
PAR_MINIAPPS =
|
||||
ifeq ($(MFEM_USE_MPI),NO)
|
||||
MINIAPPS = $(SEQ_MINIAPPS)
|
||||
else
|
||||
MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
|
||||
endif
|
||||
|
||||
.SUFFIXES:
|
||||
.SUFFIXES: .o .cpp .mk
|
||||
.PHONY: all clean clean-build clean-exec
|
||||
|
||||
# Remove built-in rules
|
||||
%: %.cpp
|
||||
%.o: %.cpp
|
||||
|
||||
# Rules for building the sources
|
||||
%.o: %.cpp $(CONFIG_MK)
|
||||
$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
|
||||
|
||||
# Replace the default implicit rule for *.cpp files
|
||||
#%: %.cpp $(UTILS)%.o $(MFEM_LIB_FILE) $(CONFIG_MK)
|
||||
# $(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
|
||||
|
||||
# Shit, I'm doing it by hand
|
||||
hptest: hptest.o exact.o util.o error.o $(MFEM_LIB_FILE)
|
||||
$(MFEM_CXX) $(MFEM_FLAGS) $(<) -o $(@) exact.o util.o error.o $(MFEM_LIBS)
|
||||
|
||||
all: $(MINIAPPS)
|
||||
|
||||
# Generate an error message if the MFEM library is not built and exit
|
||||
$(MFEM_LIB_FILE):
|
||||
$(error The MFEM library is not built)
|
||||
|
||||
clean: clean-build clean-exec
|
||||
|
||||
clean-build:
|
||||
rm -f *.o *~ $(SEQ_MINIAPPS)
|
||||
rm -rf *.dSYM *.TVD.*breakpoints
|
||||
|
||||
clean-exec:
|
||||
@rm -f partitioning.txt shaper.mesh extruder.mesh
|
||||
@rm -f optimized* perturbed*
|
||||
@@ -0,0 +1,16 @@
|
||||
set terminal postscript eps enhanced color
|
||||
set title "L-shaped domain problem"
|
||||
set size 0.9, 0.9
|
||||
set logscale y
|
||||
set format y "%1.0e"
|
||||
set grid
|
||||
set xlabel 'Square root of the number of unknowns'
|
||||
set ylabel 'H10 error'
|
||||
|
||||
plot "lshape_h1.err" using ($1**(1.0/2.0)):2 w lp pt 7 ps 0.7 linecolor 3 linetype 1 linewidth 1 title 'h-FEM, p = 1', \
|
||||
"lshape_h2.err" using ($1**(1.0/2.0)):2 w lp pt 7 ps 0.7 linecolor 2 linetype 1 linewidth 1 title 'h-FEM, p = 2', \
|
||||
"lshape_h3.err" using ($1**(1.0/2.0)):2 w lp pt 7 ps 0.7 linecolor 4 linetype 1 linewidth 1 title 'h-FEM, p = 3', \
|
||||
"../conv.err" using ($1**(1.0/2.0)):2 w lp pt 7 ps 0.7 linecolor 7 linetype 1 linewidth 1 title 'full hp-FEM', \
|
||||
|
||||
|
||||
#"lshape_hp_x_82.err" using ($1**(1.0/2.0)):2 w lp pt 7 ps 0.7 linecolor 7 linetype 1 linewidth 1 title 'hp - relaxed, 82 candidates', \
|
||||
@@ -0,0 +1,23 @@
|
||||
21 0.209401
|
||||
28 0.148039
|
||||
35 0.115693
|
||||
42 0.100031
|
||||
59 0.0791921
|
||||
92 0.0618671
|
||||
125 0.0502173
|
||||
176 0.0410161
|
||||
240 0.0341291
|
||||
381 0.0264525
|
||||
548 0.0216213
|
||||
812 0.0177322
|
||||
1134 0.0147686
|
||||
1569 0.0124069
|
||||
2206 0.0104253
|
||||
3303 0.00847206
|
||||
4641 0.00710201
|
||||
6470 0.00598583
|
||||
9259 0.00502175
|
||||
13669 0.00411453
|
||||
19226 0.00344798
|
||||
26633 0.00292862
|
||||
38436 0.00244983
|
||||
@@ -0,0 +1,32 @@
|
||||
65 0.0946127
|
||||
97 0.0600172
|
||||
129 0.0385501
|
||||
161 0.0254369
|
||||
193 0.017723
|
||||
225 0.0134895
|
||||
257 0.0113807
|
||||
369 0.00780005
|
||||
469 0.00580759
|
||||
617 0.00410815
|
||||
753 0.00322766
|
||||
981 0.00241444
|
||||
1185 0.00197134
|
||||
1661 0.00139219
|
||||
2197 0.00101529
|
||||
2849 0.000766396
|
||||
3489 0.000627526
|
||||
4377 0.000500841
|
||||
5569 0.00039321
|
||||
7449 0.00028961
|
||||
9153 0.00023255
|
||||
11453 0.000184209
|
||||
14293 0.000148004
|
||||
18361 0.000115966
|
||||
23681 9.03235e-05
|
||||
30513 6.9733e-05
|
||||
37761 5.62786e-05
|
||||
#46375 4.82672e-05
|
||||
#55563 4.28951e-05
|
||||
#65247 3.8986e-05
|
||||
#69431 3.49826e-05
|
||||
#76157 3.46286e-05
|
||||
@@ -0,0 +1,37 @@
|
||||
133 0.0556297
|
||||
208 0.0350479
|
||||
283 0.0221011
|
||||
358 0.0139625
|
||||
433 0.00885976
|
||||
508 0.0056818
|
||||
583 0.00373419
|
||||
658 0.00258199
|
||||
733 0.00194387
|
||||
808 0.00162251
|
||||
997 0.00111613
|
||||
1174 0.000833648
|
||||
1444 0.00056476
|
||||
1768 0.000380661
|
||||
2080 0.00027658
|
||||
2440 0.000206359
|
||||
2806 0.000162391
|
||||
3430 0.000116217
|
||||
4093 8.7209e-05
|
||||
4822 6.7613e-05
|
||||
5749 4.96211e-05
|
||||
6628 3.86199e-05
|
||||
8002 2.79871e-05
|
||||
9640 2.07978e-05
|
||||
11029 1.67636e-05
|
||||
13141 1.26642e-05
|
||||
15739 9.54552e-06
|
||||
18787 7.19701e-06
|
||||
22618 5.33068e-06
|
||||
26437 4.15272e-06
|
||||
31342 3.175e-06
|
||||
36847 2.46979e-06
|
||||
44701 1.84061e-06
|
||||
52903 1.429e-06
|
||||
63865 1.07194e-06
|
||||
75508 8.24746e-07
|
||||
90325 6.21383e-07
|
||||
@@ -0,0 +1,45 @@
|
||||
21 0.209401
|
||||
28 0.117024
|
||||
41 0.091868
|
||||
66 0.0621787
|
||||
107 0.0387785
|
||||
124 0.0246257
|
||||
150 0.0235255
|
||||
156 0.0176012
|
||||
182 0.0169917
|
||||
188 0.0138598
|
||||
214 0.0135543
|
||||
220 0.012057
|
||||
261 0.0100319
|
||||
291 0.00756893
|
||||
345 0.00602375
|
||||
351 0.0054997
|
||||
426 0.00394038
|
||||
432 0.00362355
|
||||
501 0.00270089
|
||||
507 0.00251869
|
||||
606 0.00176015
|
||||
612 0.00164946
|
||||
707 0.00119608
|
||||
713 0.00113174
|
||||
833 0.00078529
|
||||
839 0.000746486
|
||||
973 0.000518702
|
||||
979 0.000495442
|
||||
1127 0.000345422
|
||||
1153 0.000318529
|
||||
1299 0.000226724
|
||||
1323 0.000210516
|
||||
1469 0.000156152
|
||||
1531 0.000133411
|
||||
1688 9.78235e-05
|
||||
1730 8.95626e-05
|
||||
1953 6.22359e-05
|
||||
1995 5.70876e-05
|
||||
2210 4.01645e-05
|
||||
2252 3.70054e-05
|
||||
2491 2.55871e-05
|
||||
2533 2.3621e-05
|
||||
2770 1.65529e-05
|
||||
2812 1.53494e-05
|
||||
3075 1.08284e-05
|
||||
@@ -0,0 +1,100 @@
|
||||
#include "mfem.hpp"
|
||||
#include "util.hpp"
|
||||
|
||||
namespace mfem
|
||||
{
|
||||
|
||||
GridFunction* ProlongToMaxOrder(const GridFunction *x)
|
||||
{
|
||||
const FiniteElementSpace *fespace = x->FESpace();
|
||||
Mesh *mesh = fespace->GetMesh();
|
||||
const FiniteElementCollection *fec = fespace->FEColl();
|
||||
|
||||
// find the max order in the space
|
||||
int max_order = 1;
|
||||
for (int i = 0; i < mesh->GetNE(); i++)
|
||||
{
|
||||
max_order = std::max(fespace->GetElementOrder(i), max_order);
|
||||
}
|
||||
|
||||
// create a visualization space of max order for all elements
|
||||
FiniteElementCollection *l2fec =
|
||||
new L2_FECollection(max_order, mesh->Dimension(), BasisType::GaussLobatto);
|
||||
FiniteElementSpace *l2space = new FiniteElementSpace(mesh, l2fec);
|
||||
|
||||
IsoparametricTransformation T;
|
||||
DenseMatrix I;
|
||||
|
||||
GridFunction *prolonged_x = new GridFunction(l2space);
|
||||
|
||||
// interpolate solution vector in the larger space
|
||||
for (int i = 0; i < mesh->GetNE(); i++)
|
||||
{
|
||||
Geometry::Type geom = mesh->GetElementGeometry(i);
|
||||
T.SetIdentityTransformation(geom);
|
||||
|
||||
Array<int> dofs;
|
||||
fespace->GetElementDofs(i, dofs);
|
||||
Vector elemvect, l2vect;
|
||||
x->GetSubVector(dofs, elemvect);
|
||||
|
||||
const auto *fe = fec->GetFE(geom, fespace->GetElementOrder(i));
|
||||
const auto *l2fe = l2fec->GetFE(geom, max_order);
|
||||
|
||||
l2fe->GetTransferMatrix(*fe, T, I);
|
||||
l2space->GetElementDofs(i, dofs);
|
||||
l2vect.SetSize(dofs.Size());
|
||||
|
||||
I.Mult(elemvect, l2vect);
|
||||
prolonged_x->SetSubVector(dofs, l2vect);
|
||||
}
|
||||
|
||||
prolonged_x->MakeOwner(l2fec);
|
||||
return prolonged_x;
|
||||
}
|
||||
|
||||
|
||||
const char vishost[] = "localhost";
|
||||
const int visport = 19916;
|
||||
|
||||
void VisualizeField(socketstream &sock, GridFunction &gf, const char *title,
|
||||
const char *keys, int w, int h, int x, int y, bool vec)
|
||||
{
|
||||
Mesh &mesh = *gf.FESpace()->GetMesh();
|
||||
|
||||
bool newly_opened = false;
|
||||
int connection_failed;
|
||||
|
||||
do
|
||||
{
|
||||
if (!sock.is_open() || !sock)
|
||||
{
|
||||
sock.open(vishost, visport);
|
||||
sock.precision(8);
|
||||
newly_opened = true;
|
||||
}
|
||||
sock << "solution\n";
|
||||
|
||||
mesh.Print(sock);
|
||||
gf.Save(sock);
|
||||
|
||||
if (newly_opened)
|
||||
{
|
||||
sock << "window_title '" << title << "'\n"
|
||||
<< "window_geometry "
|
||||
<< x << " " << y << " " << w << " " << h << "\n";
|
||||
|
||||
if (keys) { sock << "keys " << keys << "\n"; }
|
||||
else { sock << "keys mAc\n"; }
|
||||
|
||||
if (vec) { sock << "vvv"; }
|
||||
sock << std::endl;
|
||||
}
|
||||
|
||||
connection_failed = !sock && !newly_opened;
|
||||
}
|
||||
while (connection_failed);
|
||||
}
|
||||
|
||||
|
||||
} // namespace mfem
|
||||
@@ -0,0 +1,20 @@
|
||||
#ifndef HPTEST_UTIL_HPP
|
||||
#define HPTEST_UTIL_HPP
|
||||
|
||||
namespace mfem
|
||||
{
|
||||
|
||||
/** Find the maximum element order of a variable-order solution 'x' and
|
||||
* prolong it to an L2 GridFunction of constant order.
|
||||
*/
|
||||
GridFunction* ProlongToMaxOrder(const GridFunction *x);
|
||||
|
||||
|
||||
void VisualizeField(socketstream &sock, GridFunction &gf, const char *title,
|
||||
const char * keys = NULL, int w = 400, int h = 400,
|
||||
int x = 0, int y = 0, bool vec = false);
|
||||
|
||||
|
||||
} // namespace mfem
|
||||
|
||||
#endif // HPTEST_UTIL_HPP
|
||||
Reference in New Issue
Block a user