1612 lines
46 KiB
C++
1612 lines
46 KiB
C++
// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
|
|
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
|
|
// reserved. See file COPYRIGHT for details.
|
|
//
|
|
// This file is part of the MFEM library. For more information and source code
|
|
// availability see http://mfem.googlecode.com.
|
|
//
|
|
// MFEM is free software; you can redistribute it and/or modify it under the
|
|
// terms of the GNU Lesser General Public License (as published by the Free
|
|
// Software Foundation) version 2.1 dated February 1999.
|
|
|
|
// Implementation of GridFunction
|
|
|
|
#include "fem.hpp"
|
|
#include <math.h>
|
|
#include <string.h>
|
|
|
|
GridFunction::GridFunction (Mesh *m, istream &input)
|
|
: Vector()
|
|
{
|
|
const int bufflen = 256;
|
|
char buff[bufflen];
|
|
int vdim;
|
|
|
|
input.getline (buff, bufflen); // 'FiniteElementSpace'
|
|
if (strcmp(buff, "FiniteElementSpace"))
|
|
mfem_error ("GridFunction::GridFunction ():"
|
|
" input stream is not a GridFunction!");
|
|
input.getline (buff, bufflen, ' '); // 'FiniteElementCollection:'
|
|
input >> ws;
|
|
input.getline (buff, bufflen);
|
|
if (!strcmp(buff, "Linear"))
|
|
fec = new LinearFECollection;
|
|
else if (!strcmp(buff, "Quadratic"))
|
|
fec = new QuadraticFECollection;
|
|
else if (!strcmp(buff, "QuadraticPos"))
|
|
fec = new QuadraticPosFECollection;
|
|
else if (!strcmp(buff, "Cubic"))
|
|
fec = new CubicFECollection;
|
|
else if (!strcmp(buff, "Const3D"))
|
|
fec = new Const3DFECollection;
|
|
else if (!strcmp(buff, "Const2D"))
|
|
fec = new Const2DFECollection;
|
|
else if (!strcmp(buff, "LinearDiscont2D"))
|
|
fec = new LinearDiscont2DFECollection;
|
|
else if (!strcmp(buff, "GaussLinearDiscont2D"))
|
|
fec = new GaussLinearDiscont2DFECollection;
|
|
else if (!strcmp(buff, "P1OnQuad"))
|
|
fec = new P1OnQuadFECollection;
|
|
else if (!strcmp(buff, "QuadraticDiscont2D"))
|
|
fec = new QuadraticDiscont2DFECollection;
|
|
else if (!strcmp(buff, "QuadraticPosDiscont2D"))
|
|
fec = new QuadraticPosDiscont2DFECollection;
|
|
else if (!strcmp(buff, "GaussQuadraticDiscont2D"))
|
|
fec = new GaussQuadraticDiscont2DFECollection;
|
|
else if (!strcmp(buff, "CubicDiscont2D"))
|
|
fec = new CubicDiscont2DFECollection;
|
|
else if (!strcmp(buff, "LinearDiscont3D"))
|
|
fec = new LinearDiscont3DFECollection;
|
|
else if (!strcmp(buff, "QuadraticDiscont3D"))
|
|
fec = new QuadraticDiscont3DFECollection;
|
|
else if (!strcmp(buff, "LinearNonConf3D"))
|
|
fec = new LinearNonConf3DFECollection;
|
|
else if (!strcmp(buff, "CrouzeixRaviart"))
|
|
fec = new CrouzeixRaviartFECollection;
|
|
else if (!strcmp(buff, "ND1_3D"))
|
|
fec = new ND1_3DFECollection;
|
|
else
|
|
mfem_error ("GridFunction::GridFunction (): "
|
|
"Unknown FiniteElementCollection!");
|
|
input.getline (buff, bufflen, ' '); // 'VDim:'
|
|
input >> vdim;
|
|
input.getline (buff, bufflen, ' '); // 'Ordering:'
|
|
int ordering;
|
|
input >> ordering;
|
|
input.getline (buff, bufflen); // read the empty line
|
|
fes = new FiniteElementSpace (m, fec, vdim, ordering);
|
|
Vector::Load (input, fes -> GetVSize());
|
|
}
|
|
|
|
GridFunction::~GridFunction ()
|
|
{
|
|
if (fec)
|
|
{
|
|
delete fes;
|
|
delete fec;
|
|
}
|
|
}
|
|
|
|
void GridFunction::Update(FiniteElementSpace *f)
|
|
{
|
|
if (fec)
|
|
{
|
|
delete fes;
|
|
delete fec;
|
|
fec = NULL;
|
|
}
|
|
fes = f;
|
|
SetSize(fes->GetVSize());
|
|
}
|
|
|
|
|
|
void GridFunction::Update(FiniteElementSpace *f, Vector &v, int v_offset)
|
|
{
|
|
if (fec)
|
|
{
|
|
delete fes;
|
|
delete fec;
|
|
fec = NULL;
|
|
}
|
|
fes = f;
|
|
SetDataAndSize((double *)v + v_offset, fes->GetVSize());
|
|
}
|
|
|
|
int GridFunction::VectorDim() const
|
|
{
|
|
const FiniteElement *fe = fes -> GetFE(0);
|
|
|
|
if (fe -> GetRangeType() == FiniteElement::SCALAR)
|
|
return fes -> GetVDim();
|
|
return fe -> GetDim();
|
|
}
|
|
|
|
void GridFunction::GetNodalValues (int i, Array<double> & nval, int vdim) const
|
|
{
|
|
Array<int> dofs;
|
|
|
|
int k;
|
|
|
|
fes -> GetElementVDofs (i, dofs);
|
|
const FiniteElement * FElem = fes -> GetFE (i);
|
|
const IntegrationRule * ElemVert =
|
|
Geometries.GetVertices (FElem -> GetGeomType ());
|
|
int dof = FElem -> GetDof ();
|
|
Vector DofVal (dof);
|
|
int n = ElemVert -> GetNPoints ();
|
|
nval.SetSize (n);
|
|
vdim--;
|
|
for (k = 0; k < n; k++)
|
|
{
|
|
FElem -> CalcShape (ElemVert -> IntPoint (k), DofVal);
|
|
nval[k] = 0.0;
|
|
for (int j = 0; j < dof; j++)
|
|
if (dofs[dof*vdim+j] >= 0)
|
|
nval[k] += DofVal (j) * data[dofs[dof*vdim+j]];
|
|
else
|
|
nval[k] -= DofVal (j) * data[-1-dofs[dof*vdim+j]];
|
|
}
|
|
}
|
|
|
|
|
|
double GridFunction::GetValue (int i, const IntegrationPoint &ip, int vdim)
|
|
const
|
|
{
|
|
Array<int> dofs;
|
|
fes -> GetElementDofs (i, dofs);
|
|
fes -> DofsToVDofs (vdim-1, dofs);
|
|
Vector DofVal (dofs.Size()), LocVec;
|
|
fes -> GetFE (i) -> CalcShape (ip, DofVal);
|
|
GetSubVector (dofs, LocVec);
|
|
|
|
return (DofVal * LocVec);
|
|
}
|
|
|
|
void GridFunction::GetVectorValue (int i, const IntegrationPoint &ip,
|
|
Vector &val) const
|
|
{
|
|
const FiniteElement *FElem = fes -> GetFE (i);
|
|
int dof = FElem -> GetDof();
|
|
Array<int> vdofs;
|
|
fes -> GetElementVDofs (i, vdofs);
|
|
Vector loc_data;
|
|
GetSubVector (vdofs, loc_data);
|
|
if (FElem -> GetRangeType() == FiniteElement::SCALAR)
|
|
{
|
|
Vector shape (dof);
|
|
FElem -> CalcShape (ip, shape);
|
|
int vdim = fes -> GetVDim();
|
|
val.SetSize (vdim);
|
|
for (int k = 0; k < vdim; k++)
|
|
{
|
|
val(k) = shape * ((const double *)loc_data + dof * k);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
int dim = FElem -> GetDim();
|
|
DenseMatrix vshape (dof, dim);
|
|
ElementTransformation *Tr = fes -> GetElementTransformation (i);
|
|
Tr -> SetIntPoint (&ip);
|
|
FElem -> CalcVShape (*Tr, vshape);
|
|
val.SetSize (dim);
|
|
vshape.MultTranspose (loc_data, val);
|
|
}
|
|
}
|
|
|
|
void GridFunction::GetValues (int i, const IntegrationRule &ir, Vector &vals,
|
|
DenseMatrix &tr, int vdim)
|
|
const
|
|
{
|
|
Array<int> dofs;
|
|
|
|
int k, n;
|
|
|
|
n = ir.GetNPoints ();
|
|
vals.SetSize (n);
|
|
fes -> GetElementVDofs (i, dofs);
|
|
const FiniteElement * FElem = fes -> GetFE (i);
|
|
ElementTransformation *ET;
|
|
ET = fes -> GetElementTransformation (i);
|
|
ET -> Transform (ir, tr);
|
|
int dof = FElem -> GetDof ();
|
|
Vector DofVal (dof);
|
|
vdim--;
|
|
for (k = 0; k < n; k++)
|
|
{
|
|
FElem -> CalcShape (ir.IntPoint (k), DofVal);
|
|
vals(k) = 0.0;
|
|
for (int j = 0; j < dof; j++)
|
|
if (dofs[dof*vdim+j] >= 0)
|
|
vals(k) += DofVal (j) * data[dofs[dof*vdim+j]];
|
|
else
|
|
vals(k) -= DofVal (j) * data[-1-dofs[dof*vdim+j]];
|
|
}
|
|
}
|
|
|
|
int GridFunction::GetFaceValues (int i, int side, const IntegrationRule &ir,
|
|
Vector &vals, DenseMatrix &tr,
|
|
int vdim) const
|
|
{
|
|
int n, di;
|
|
FaceElementTransformations *Transf;
|
|
|
|
n = ir.GetNPoints();
|
|
IntegrationRule eir (n); // ---
|
|
Transf = fes -> GetMesh() -> GetFaceElementTransformations (i);
|
|
if (Transf -> Elem2No < 0)
|
|
di = 0;
|
|
else
|
|
if ( fes -> GetAttribute (Transf -> Elem1No) <=
|
|
fes -> GetAttribute (Transf -> Elem2No) )
|
|
di = 0;
|
|
else
|
|
di = 1;
|
|
if (side == 2)
|
|
side = di;
|
|
if (side == 0)
|
|
{
|
|
Transf -> Loc1.Transform (ir, eir);
|
|
GetValues (Transf -> Elem1No, eir, vals, tr, vdim);
|
|
}
|
|
else
|
|
{
|
|
Transf -> Loc2.Transform (ir, eir);
|
|
GetValues (Transf -> Elem2No, eir, vals, tr, vdim);
|
|
}
|
|
|
|
return di;
|
|
}
|
|
|
|
void GridFunction::GetVectorValues (int i, const IntegrationRule &ir,
|
|
DenseMatrix &vals, DenseMatrix &tr)
|
|
const
|
|
{
|
|
Array<int> dofs;
|
|
|
|
int k, n, vdim;
|
|
|
|
n = ir.GetNPoints ();
|
|
vdim = fes->GetVDim();
|
|
vals.SetSize (vdim, n);
|
|
fes -> GetElementVDofs (i, dofs);
|
|
const FiniteElement * FElem = fes -> GetFE (i);
|
|
ElementTransformation *ET;
|
|
ET = fes -> GetElementTransformation (i);
|
|
ET -> Transform (ir, tr);
|
|
int dof = FElem -> GetDof ();
|
|
Vector DofVal (dof);
|
|
for (k = 0; k < n; k++)
|
|
{
|
|
FElem -> CalcShape (ir.IntPoint (k), DofVal);
|
|
for (int d = 0; d < vdim; d++)
|
|
{
|
|
double v = 0.0;
|
|
for (int j = 0; j < dof; j++)
|
|
{
|
|
int ind = dofs[dof*d+j];
|
|
if (ind >= 0)
|
|
v += DofVal (j) * data[ind];
|
|
else
|
|
v -= DofVal (j) * data[-1-ind];
|
|
}
|
|
vals(d,k) = v;
|
|
}
|
|
}
|
|
}
|
|
|
|
int GridFunction::GetFaceVectorValues (
|
|
int i, int side, const IntegrationRule &ir,
|
|
DenseMatrix &vals, DenseMatrix &tr) const
|
|
{
|
|
int n, di;
|
|
FaceElementTransformations *Transf;
|
|
|
|
n = ir.GetNPoints();
|
|
IntegrationRule eir (n); // ---
|
|
Transf = fes -> GetMesh() -> GetFaceElementTransformations (i);
|
|
if (Transf -> Elem2No < 0)
|
|
di = 0;
|
|
else
|
|
if ( fes -> GetAttribute (Transf -> Elem1No) <=
|
|
fes -> GetAttribute (Transf -> Elem2No) )
|
|
di = 0;
|
|
else
|
|
di = 1;
|
|
if (side == 2)
|
|
side = di;
|
|
if (side == 0)
|
|
{
|
|
Transf -> Loc1.Transform (ir, eir);
|
|
GetVectorValues (Transf -> Elem1No, eir, vals, tr);
|
|
}
|
|
else
|
|
{
|
|
Transf -> Loc2.Transform (ir, eir);
|
|
GetVectorValues (Transf -> Elem2No, eir, vals, tr);
|
|
}
|
|
|
|
return di;
|
|
}
|
|
|
|
void GridFunction::GetValuesFrom (GridFunction &orig_func)
|
|
{
|
|
// Without averaging ...
|
|
|
|
FiniteElementSpace *orig_fes = orig_func.FESpace();
|
|
Array<int> vdofs, orig_vdofs;
|
|
Vector shape, loc_values, orig_loc_values;
|
|
int i, j, d, ne, dof, odof, vdim;
|
|
|
|
ne = fes -> GetNE();
|
|
vdim = fes -> GetVDim();
|
|
for (i = 0; i < ne; i++)
|
|
{
|
|
fes -> GetElementVDofs (i, vdofs);
|
|
orig_fes -> GetElementVDofs (i, orig_vdofs);
|
|
orig_func.GetSubVector(orig_vdofs, orig_loc_values);
|
|
const FiniteElement *fe = fes -> GetFE (i);
|
|
const FiniteElement *orig_fe = orig_fes -> GetFE (i);
|
|
dof = fe -> GetDof();
|
|
odof = orig_fe -> GetDof();
|
|
loc_values.SetSize (dof * vdim);
|
|
shape.SetSize (odof);
|
|
const IntegrationRule &ir = fe -> GetNodes();
|
|
for (j = 0; j < dof; j++)
|
|
{
|
|
const IntegrationPoint &ip = ir.IntPoint (j);
|
|
orig_fe -> CalcShape (ip, shape);
|
|
for (d = 0; d < vdim; d++)
|
|
{
|
|
loc_values(d*dof+j) =
|
|
shape * ((const double *)orig_loc_values + d * odof) ;
|
|
}
|
|
}
|
|
SetSubVector (vdofs, loc_values);
|
|
}
|
|
}
|
|
|
|
void GridFunction::GetBdrValuesFrom (GridFunction &orig_func)
|
|
{
|
|
// Without averaging ...
|
|
|
|
FiniteElementSpace *orig_fes = orig_func.FESpace();
|
|
Array<int> vdofs, orig_vdofs;
|
|
Vector shape, loc_values, orig_loc_values;
|
|
int i, j, d, nbe, dof, odof, vdim;
|
|
|
|
nbe = fes -> GetNBE();
|
|
vdim = fes -> GetVDim();
|
|
for (i = 0; i < nbe; i++)
|
|
{
|
|
fes -> GetBdrElementVDofs (i, vdofs);
|
|
orig_fes -> GetBdrElementVDofs (i, orig_vdofs);
|
|
orig_func.GetSubVector(orig_vdofs, orig_loc_values);
|
|
const FiniteElement *fe = fes -> GetBE (i);
|
|
const FiniteElement *orig_fe = orig_fes -> GetBE (i);
|
|
dof = fe -> GetDof();
|
|
odof = orig_fe -> GetDof();
|
|
loc_values.SetSize (dof * vdim);
|
|
shape.SetSize (odof);
|
|
const IntegrationRule &ir = fe -> GetNodes();
|
|
for (j = 0; j < dof; j++)
|
|
{
|
|
const IntegrationPoint &ip = ir.IntPoint (j);
|
|
orig_fe -> CalcShape (ip, shape);
|
|
for (d = 0; d < vdim; d++)
|
|
{
|
|
loc_values(d*dof+j) =
|
|
shape * ((const double *)orig_loc_values + d * odof);
|
|
}
|
|
}
|
|
SetSubVector (vdofs, loc_values);
|
|
}
|
|
}
|
|
|
|
void GridFunction::GetVectorFieldValues (
|
|
int i, const IntegrationRule &ir, DenseMatrix &vals,
|
|
DenseMatrix &tr, int comp) const
|
|
{
|
|
Array<int> vdofs;
|
|
ElementTransformation *transf;
|
|
|
|
int d, j, k, n, dim, dof, ind;
|
|
|
|
n = ir.GetNPoints ();
|
|
fes -> GetElementVDofs (i, vdofs);
|
|
const FiniteElement *fe = fes -> GetFE(i);
|
|
dof = fe -> GetDof ();
|
|
dim = fe -> GetDim();
|
|
int *dofs = &vdofs[comp*dof];
|
|
transf = fes -> GetElementTransformation (i);
|
|
transf -> Transform (ir, tr);
|
|
vals.SetSize (n, dim);
|
|
DenseMatrix vshape (dof, dim);
|
|
double a;
|
|
for (k = 0; k < n; k++)
|
|
{
|
|
const IntegrationPoint &ip = ir.IntPoint (k);
|
|
transf -> SetIntPoint (&ip);
|
|
fe -> CalcVShape (*transf, vshape);
|
|
for (d = 0; d < dim; d++)
|
|
{
|
|
a = 0.0;
|
|
for (j = 0; j < dof; j++)
|
|
if ( (ind=dofs[j]) >= 0 )
|
|
a += vshape(j, d) * data[ind];
|
|
else
|
|
a -= vshape(j, d) * data[-1-ind];
|
|
vals(k, d) = a;
|
|
}
|
|
}
|
|
}
|
|
|
|
void GridFunction::ReorderByNodes()
|
|
{
|
|
if (fes->GetOrdering() == Ordering::byNODES)
|
|
return;
|
|
|
|
int i, j, k;
|
|
int vdim = fes -> GetVDim();
|
|
int ndofs = fes -> GetNDofs();
|
|
double *temp = new double[size];
|
|
|
|
k = 0;
|
|
for (j = 0; j < ndofs; j++)
|
|
for (i = 0; i < vdim; i++)
|
|
temp[j+i*ndofs] = data[k++];
|
|
|
|
for (i = 0; i < size; i++)
|
|
data[i] = temp[i];
|
|
|
|
delete [] temp;
|
|
}
|
|
|
|
void GridFunction::GetVectorFieldNodalValues (Vector &val, int comp) const
|
|
{
|
|
int i, k;
|
|
Array<int> overlap (fes -> GetNV());
|
|
Array<int> vertices;
|
|
DenseMatrix vals, tr;
|
|
|
|
val.SetSize (overlap.Size());
|
|
for (i = 0; i < overlap.Size(); i++)
|
|
overlap[i] = 0;
|
|
val = 0.0;
|
|
|
|
comp--;
|
|
for (i = 0; i < fes -> GetNE(); i++)
|
|
{
|
|
const IntegrationRule *ir =
|
|
Geometries.GetVertices (fes -> GetFE(i) -> GetGeomType ());
|
|
fes -> GetElementVertices (i, vertices);
|
|
GetVectorFieldValues (i, *ir, vals, tr);
|
|
for (k = 0; k < ir -> GetNPoints(); k++)
|
|
{
|
|
val(vertices[k]) += vals(k, comp);
|
|
overlap[vertices[k]]++;
|
|
}
|
|
}
|
|
|
|
for (i = 0; i < overlap.Size(); i++)
|
|
val(i) /= overlap[i];
|
|
}
|
|
|
|
void GridFunction::ProjectVectorFieldOn (GridFunction &vec_field, int comp)
|
|
{
|
|
FiniteElementSpace * new_fes = vec_field.FESpace();
|
|
|
|
int d, i, k, ind, dof;
|
|
Array<int> overlap (new_fes -> GetVSize());
|
|
Array<int> new_vdofs;
|
|
DenseMatrix vals, tr;
|
|
|
|
for (i = 0; i < overlap.Size(); i++)
|
|
overlap[i] = 0;
|
|
vec_field = 0.0;
|
|
|
|
for (i = 0; i < new_fes -> GetNE(); i++)
|
|
{
|
|
const FiniteElement *fe = new_fes -> GetFE(i);
|
|
const IntegrationRule &ir = fe -> GetNodes();
|
|
GetVectorFieldValues (i, ir, vals, tr, comp);
|
|
new_fes -> GetElementVDofs (i, new_vdofs);
|
|
dof = fe -> GetDof();
|
|
for (d = 0; d < fe -> GetDim(); d++)
|
|
for (k = 0; k < dof; k++)
|
|
{
|
|
if ( (ind=new_vdofs[dof*d+k]) < 0 )
|
|
ind = -1-ind, vals(k, d) = - vals(k, d);
|
|
vec_field(ind) += vals(k, d);
|
|
overlap[ind]++;
|
|
}
|
|
}
|
|
|
|
for (i = 0; i < overlap.Size(); i++)
|
|
vec_field(i) /= overlap[i];
|
|
}
|
|
|
|
void GridFunction::GetDerivative(int comp, int der_comp, GridFunction &der)
|
|
{
|
|
FiniteElementSpace * der_fes = der.FESpace();
|
|
ElementTransformation * transf;
|
|
Array<int> overlap (der_fes -> GetVSize());
|
|
Array<int> der_dofs, vdofs;
|
|
DenseMatrix dshape, inv_jac;
|
|
Vector pt_grad, loc_func;
|
|
int i, j, k, dim, dof, der_dof, ind;
|
|
double a;
|
|
|
|
for (i = 0; i < overlap.Size(); i++)
|
|
overlap[i] = 0;
|
|
der = 0.0;
|
|
|
|
comp--;
|
|
for (i = 0; i < der_fes -> GetNE(); i++)
|
|
{
|
|
const FiniteElement *der_fe = der_fes -> GetFE(i);
|
|
const FiniteElement *fe = fes -> GetFE(i);
|
|
const IntegrationRule &ir = der_fe -> GetNodes();
|
|
der_fes -> GetElementDofs (i, der_dofs);
|
|
fes -> GetElementVDofs (i, vdofs);
|
|
dim = fe -> GetDim();
|
|
dof = fe -> GetDof();
|
|
der_dof = der_fe -> GetDof();
|
|
dshape.SetSize (dof, dim);
|
|
inv_jac.SetSize (dim);
|
|
pt_grad.SetSize (dim);
|
|
loc_func.SetSize (dof);
|
|
transf = fes -> GetElementTransformation (i);
|
|
for (j = 0; j < dof; j++)
|
|
loc_func(j) = ( (ind=vdofs[comp*dof+j]) >= 0 ) ?
|
|
(data[ind]) : (-data[-1-ind]);
|
|
for (k = 0; k < der_dof; k++)
|
|
{
|
|
const IntegrationPoint &ip = ir.IntPoint (k);
|
|
fe -> CalcDShape (ip, dshape);
|
|
dshape.MultTranspose (loc_func, pt_grad);
|
|
transf -> SetIntPoint (&ip);
|
|
CalcInverse (transf -> Jacobian(), inv_jac);
|
|
a = 0.0;
|
|
for (j = 0; j < dim; j++)
|
|
a += inv_jac(j, der_comp) * pt_grad(j);
|
|
der(der_dofs[k]) += a;
|
|
overlap[der_dofs[k]]++;
|
|
}
|
|
}
|
|
|
|
for (i = 0; i < overlap.Size(); i++)
|
|
der(i) /= overlap[i];
|
|
}
|
|
|
|
|
|
void GridFunction::GetVectorGradientHat(
|
|
ElementTransformation &T, DenseMatrix &gh)
|
|
{
|
|
int elNo = T.ElementNo;
|
|
const FiniteElement *FElem = fes -> GetFE (elNo);
|
|
int dim = FElem -> GetDim(), dof = FElem -> GetDof();
|
|
Array<int> vdofs;
|
|
fes -> GetElementVDofs (elNo, vdofs);
|
|
Vector loc_data;
|
|
GetSubVector (vdofs, loc_data);
|
|
// assuming scalar FE
|
|
DenseMatrix dshape (dof, dim);
|
|
FElem -> CalcDShape (T.GetIntPoint(), dshape);
|
|
gh.SetSize(dim);
|
|
for (int i = 0; i < dim; i++)
|
|
for (int j = 0; j < dim; j++)
|
|
{
|
|
double gij = 0.0;
|
|
for (int k = 0; k < dof; k++)
|
|
gij += loc_data(i * dof + k) * dshape(k, j);
|
|
gh(i, j) = gij;
|
|
}
|
|
}
|
|
|
|
double GridFunction::GetDivergence(ElementTransformation &tr)
|
|
{
|
|
DenseMatrix grad_hat;
|
|
GetVectorGradientHat(tr, grad_hat);
|
|
int dim = grad_hat.Size();
|
|
DenseMatrix Jinv(dim);
|
|
CalcInverse(tr.Jacobian(), Jinv);
|
|
double div_v = 0.0;
|
|
for (int i = 0; i < dim; i++)
|
|
for (int j = 0; j < dim; j++)
|
|
div_v += grad_hat(i, j) * Jinv(j, i);
|
|
|
|
return div_v;
|
|
}
|
|
|
|
void GridFunction::GetGradient(ElementTransformation &tr, Vector &grad)
|
|
{
|
|
mfem_error("GridFunction::GetGradient(...) is not implemented!");
|
|
}
|
|
|
|
void GridFunction::GetGradients(const int elem, const IntegrationRule &ir,
|
|
DenseMatrix &grad)
|
|
{
|
|
const FiniteElement *fe = fes->GetFE(elem);
|
|
ElementTransformation *Tr = fes->GetElementTransformation(elem);
|
|
DenseMatrix dshape(fe->GetDof(), fe->GetDim());
|
|
DenseMatrix Jinv(fe->GetDim());
|
|
Vector lval, gh(fe->GetDim()), gcol;
|
|
Array<int> dofs;
|
|
fes->GetElementDofs(elem, dofs);
|
|
GetSubVector(dofs, lval);
|
|
grad.SetSize(fe->GetDim(), ir.GetNPoints());
|
|
for (int i = 0; i < ir.GetNPoints(); i++)
|
|
{
|
|
const IntegrationPoint &ip = ir.IntPoint(i);
|
|
fe->CalcDShape(ip, dshape);
|
|
dshape.MultTranspose(lval, gh);
|
|
Tr->SetIntPoint(&ip);
|
|
grad.GetColumnReference(i, gcol);
|
|
CalcInverse(Tr->Jacobian(), Jinv);
|
|
Jinv.MultTranspose(gh, gcol);
|
|
}
|
|
}
|
|
|
|
void GridFunction::GetVectorGradient(
|
|
ElementTransformation &tr, DenseMatrix &grad)
|
|
{
|
|
DenseMatrix grad_hat;
|
|
GetVectorGradientHat(tr, grad_hat);
|
|
DenseMatrix Jinv(grad_hat.Size());
|
|
CalcInverse(tr.Jacobian(), Jinv);
|
|
grad.SetSize(grad_hat.Size());
|
|
Mult(grad_hat, Jinv, grad);
|
|
}
|
|
|
|
void GridFunction::GetElementAverages (GridFunction &avgs)
|
|
{
|
|
MassIntegrator Mi;
|
|
DenseMatrix loc_mass;
|
|
Array<int> te_dofs, tr_dofs;
|
|
Vector loc_avgs, loc_this;
|
|
Vector int_psi (avgs.Size());
|
|
|
|
avgs = 0.0;
|
|
int_psi = 0.0;
|
|
for (int i = 0; i < fes -> GetNE(); i++)
|
|
{
|
|
Mi.AssembleElementMatrix2 (*fes -> GetFE (i),
|
|
*avgs.FESpace() -> GetFE (i),
|
|
*fes -> GetElementTransformation (i),
|
|
loc_mass);
|
|
fes -> GetElementDofs (i, tr_dofs);
|
|
avgs.FESpace() -> GetElementDofs (i, te_dofs);
|
|
GetSubVector (tr_dofs, loc_this);
|
|
loc_avgs.SetSize (te_dofs.Size());
|
|
loc_mass.Mult (loc_this, loc_avgs);
|
|
avgs.AddElementVector (te_dofs, loc_avgs);
|
|
loc_this = 1.0; // assume the local basis for 'this' sums to 1
|
|
loc_mass.Mult (loc_this, loc_avgs);
|
|
int_psi.AddElementVector (te_dofs, loc_avgs);
|
|
}
|
|
for (int i = 0; i < avgs.Size(); i++)
|
|
avgs(i) /= int_psi(i);
|
|
}
|
|
|
|
void GridFunction::GetNodalValues (Vector & nval, int vdim) const
|
|
{
|
|
int i, j;
|
|
Array<int> vertices;
|
|
Array<double> values;
|
|
Array<int> overlap(fes->GetNV());
|
|
nval.SetSize (fes -> GetNV());
|
|
|
|
for (i = 0; i < overlap.Size(); i++)
|
|
{
|
|
nval(i) = 0.0;
|
|
overlap[i] = 0;
|
|
}
|
|
for (i = 0; i < fes -> GetNE(); i++)
|
|
{
|
|
fes -> GetElementVertices (i, vertices);
|
|
GetNodalValues (i, values, vdim);
|
|
for (j = 0; j < vertices.Size(); j++)
|
|
{
|
|
nval(vertices[j]) += values[j];
|
|
overlap[vertices[j]]++;
|
|
}
|
|
}
|
|
for (i = 0; i < overlap.Size(); i++)
|
|
nval(i) /= overlap[i];
|
|
}
|
|
|
|
void GridFunction::ProjectCoefficient (Coefficient &coeff)
|
|
{
|
|
int i;
|
|
Array<int> vdofs;
|
|
Vector vals;
|
|
|
|
DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
|
|
|
|
if (delta_c == NULL)
|
|
{
|
|
for (i = 0; i < fes -> GetNE(); i++)
|
|
{
|
|
fes -> GetElementVDofs (i, vdofs);
|
|
vals.SetSize (vdofs.Size());
|
|
fes -> GetFE (i) -> Project (
|
|
coeff, *fes -> GetElementTransformation(i), vals);
|
|
SetSubVector (vdofs, vals);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
Mesh *mesh = fes->GetMesh();
|
|
const int dim = mesh->Dimension();
|
|
const double *center = delta_c->Center();
|
|
const double *vert = mesh->GetVertex(0);
|
|
double min_dist, dist;
|
|
int v_idx = 0;
|
|
|
|
// find the vertex closest to the center of the delta function
|
|
min_dist = Distance(center, vert, dim);
|
|
for (int i = 0; i < mesh->GetNV(); i++)
|
|
{
|
|
vert = mesh->GetVertex(i);
|
|
dist = Distance(center, vert, dim);
|
|
if (dist < min_dist)
|
|
{
|
|
min_dist = dist;
|
|
v_idx = i;
|
|
}
|
|
}
|
|
|
|
(*this) = 0.0;
|
|
|
|
// find the elements that have 'v_idx' as a vertex
|
|
MassIntegrator Mi(*delta_c->Weight());
|
|
DenseMatrix loc_mass;
|
|
Array<int> vertices;
|
|
Vector loc_mass_vals;
|
|
double integral = 0.0;
|
|
for (int i = 0; i < mesh->GetNE(); i++)
|
|
{
|
|
mesh->GetElementVertices(i, vertices);
|
|
for (int j = 0; j < vertices.Size(); j++)
|
|
if (vertices[j] == v_idx)
|
|
{
|
|
const FiniteElement *fe = fes->GetFE(i);
|
|
Mi.AssembleElementMatrix(*fe,
|
|
*fes->GetElementTransformation(i),
|
|
loc_mass);
|
|
vals.SetSize(fe->GetDof());
|
|
fe->ProjectDelta(j, vals);
|
|
fes->GetElementVDofs(i, vdofs);
|
|
SetSubVector(vdofs, vals);
|
|
loc_mass_vals.SetSize(vals.Size());
|
|
loc_mass.Mult(vals, loc_mass_vals);
|
|
for (int k = 0; k < loc_mass_vals.Size(); k++)
|
|
integral += loc_mass_vals(k);
|
|
break;
|
|
}
|
|
}
|
|
|
|
(*this) *= (delta_c->Scale() / integral);
|
|
}
|
|
}
|
|
|
|
void GridFunction::ProjectCoefficient (
|
|
Coefficient &coeff, Array<int> &dofs, int vd)
|
|
{
|
|
int el = -1;
|
|
ElementTransformation *T = NULL;
|
|
const FiniteElement *fe = NULL;
|
|
|
|
for (int i = 0; i < dofs.Size(); i++)
|
|
{
|
|
int dof = dofs[i], j = fes -> GetElementForDof (dof);
|
|
if (el != j)
|
|
{
|
|
el = j;
|
|
T = fes -> GetElementTransformation (el);
|
|
fe = fes -> GetFE (el);
|
|
}
|
|
int vdof = fes -> DofToVDof (dof, vd);
|
|
int ld = fes -> GetLocalDofForDof (dof);
|
|
const IntegrationPoint &ip = fe -> GetNodes().IntPoint(ld);
|
|
T->SetIntPoint(&ip);
|
|
(*this)(vdof) = coeff.Eval (*T, ip);
|
|
}
|
|
}
|
|
|
|
void GridFunction::ProjectCoefficient (VectorCoefficient &vcoeff)
|
|
{
|
|
int i;
|
|
Array<int> vdofs;
|
|
Vector vals;
|
|
|
|
for (i = 0; i < fes -> GetNE(); i++)
|
|
{
|
|
fes -> GetElementVDofs (i, vdofs);
|
|
vals.SetSize (vdofs.Size());
|
|
fes -> GetFE(i) -> Project (
|
|
vcoeff, *fes->GetElementTransformation(i), vals);
|
|
SetSubVector (vdofs, vals);
|
|
}
|
|
}
|
|
|
|
void GridFunction::ProjectCoefficient (Coefficient *coeff[])
|
|
{
|
|
int i, j, fdof, d, ind;
|
|
double val;
|
|
const FiniteElement *fe;
|
|
ElementTransformation *transf;
|
|
Array<int> vdofs;
|
|
|
|
for (i = 0; i < fes->GetNE(); i++)
|
|
{
|
|
fe = fes->GetFE(i);
|
|
fdof = fe->GetDof();
|
|
transf = fes->GetElementTransformation(i);
|
|
const IntegrationRule &ir = fe->GetNodes();
|
|
fes->GetElementVDofs(i, vdofs);
|
|
for (j = 0; j < fdof; j++)
|
|
{
|
|
const IntegrationPoint &ip = ir.IntPoint(j);
|
|
transf->SetIntPoint(&ip);
|
|
for (d = 0; d < fes->GetVDim(); d++)
|
|
{
|
|
val = coeff[d]->Eval(*transf, ip);
|
|
if ( (ind = vdofs[fdof*d+j]) < 0 )
|
|
val = -val, ind = -1-ind;
|
|
(*this)(ind) = val;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void GridFunction::ProjectBdrCoefficient (
|
|
Coefficient *coeff[], Array<int> &attr )
|
|
{
|
|
int i, j, fdof, d, ind;
|
|
double val;
|
|
const FiniteElement *fe;
|
|
ElementTransformation *transf;
|
|
Array<int> vdofs;
|
|
|
|
for (i = 0; i < fes->GetNBE(); i++)
|
|
{
|
|
if ( attr[fes->GetBdrAttribute(i)-1] )
|
|
{
|
|
fe = fes->GetBE(i);
|
|
fdof = fe->GetDof();
|
|
transf = fes->GetBdrElementTransformation(i);
|
|
const IntegrationRule &ir = fe->GetNodes();
|
|
fes->GetBdrElementVDofs(i, vdofs);
|
|
for (j = 0; j < fdof; j++)
|
|
{
|
|
const IntegrationPoint &ip = ir.IntPoint(j);
|
|
transf->SetIntPoint(&ip);
|
|
for (d = 0; d < fes->GetVDim(); d++)
|
|
{
|
|
val = coeff[d]->Eval(*transf, ip);
|
|
if ( (ind = vdofs[fdof*d+j]) < 0 )
|
|
val = -val, ind = -1-ind;
|
|
(*this)(ind) = val;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
double GridFunction::ComputeL2Error (
|
|
Coefficient *exsol[], const IntegrationRule *irs[]) const
|
|
{
|
|
double error = 0.0, a;
|
|
const FiniteElement *fe;
|
|
ElementTransformation *transf;
|
|
Vector shape;
|
|
Array<int> vdofs;
|
|
int fdof, d, i, intorder, j, k;
|
|
|
|
for (i = 0; i < fes->GetNE(); i++)
|
|
{
|
|
fe = fes->GetFE(i);
|
|
fdof = fe->GetDof();
|
|
transf = fes->GetElementTransformation(i);
|
|
shape.SetSize(fdof);
|
|
intorder = fe->GetOrder()+2; // <----------
|
|
const IntegrationRule *ir;
|
|
if (irs)
|
|
ir = irs[fe->GetGeomType()];
|
|
else
|
|
ir = &(IntRules.Get(fe->GetGeomType(), intorder));
|
|
fes->GetElementVDofs (i, vdofs);
|
|
for (j = 0; j < ir->GetNPoints(); j++)
|
|
{
|
|
const IntegrationPoint &ip = ir->IntPoint(j);
|
|
fe->CalcShape(ip, shape);
|
|
for (d = 0; d < fes->GetVDim(); d++)
|
|
{
|
|
a = 0;
|
|
for (k = 0; k < fdof; k++)
|
|
if (vdofs[fdof*d+k] >= 0)
|
|
a += (*this)(vdofs[fdof*d+k]) * shape(k);
|
|
else
|
|
a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
|
|
transf->SetIntPoint (&ip);
|
|
a -= exsol[d]->Eval(*transf, ip);
|
|
error += ip.weight * transf->Weight() * a * a;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (error < 0.0)
|
|
return -sqrt(-error);
|
|
return sqrt(error);
|
|
}
|
|
|
|
double GridFunction::ComputeL2Error (
|
|
VectorCoefficient &exsol, const IntegrationRule *irs[],
|
|
Array<int> *elems) const
|
|
{
|
|
double error = 0.0, a;
|
|
const FiniteElement *fe;
|
|
ElementTransformation *transf;
|
|
Vector shape, sol (exsol.GetVDim());
|
|
Array<int> vdofs;
|
|
int fdof, d, i, intorder, j, k;
|
|
|
|
for (i = 0; i < fes->GetNE(); i++)
|
|
{
|
|
if (elems != NULL && (*elems)[i] == 0) continue;
|
|
fe = fes->GetFE(i);
|
|
fdof = fe->GetDof();
|
|
transf = fes->GetElementTransformation(i);
|
|
shape.SetSize(fdof);
|
|
intorder = fe->GetOrder()+2; // <----------
|
|
const IntegrationRule *ir;
|
|
if (irs)
|
|
ir = irs[fe->GetGeomType()];
|
|
else
|
|
ir = &(IntRules.Get(fe->GetGeomType(), intorder));
|
|
fes->GetElementVDofs (i, vdofs);
|
|
for (j = 0; j < ir->GetNPoints(); j++)
|
|
{
|
|
const IntegrationPoint &ip = ir->IntPoint(j);
|
|
fe->CalcShape(ip, shape);
|
|
transf->SetIntPoint(&ip);
|
|
exsol.Eval (sol, *transf, ip);
|
|
for (d = 0; d < fes->GetVDim(); d++)
|
|
{
|
|
a = 0;
|
|
for (k = 0; k < fdof; k++)
|
|
if (vdofs[fdof*d+k] >= 0)
|
|
a += (*this)(vdofs[fdof*d+k]) * shape(k);
|
|
else
|
|
a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
|
|
a -= sol(d);
|
|
error += ip.weight * transf->Weight() * a * a;
|
|
}
|
|
}
|
|
}
|
|
|
|
if (error < 0.0)
|
|
return -sqrt(-error);
|
|
return sqrt(error);
|
|
}
|
|
|
|
double GridFunction::ComputeH1Error (
|
|
Coefficient *exsol, VectorCoefficient *exgrad,
|
|
Coefficient *ell_coeff, double Nu, int norm_type) const
|
|
{
|
|
// assuming vdim is 1
|
|
int i, fdof, dim, intorder, j, k;
|
|
Mesh *mesh;
|
|
const FiniteElement *fe;
|
|
ElementTransformation *transf;
|
|
FaceElementTransformations *face_elem_transf;
|
|
Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
|
|
DenseMatrix dshape, dshapet, Jinv;
|
|
Array<int> vdofs;
|
|
IntegrationPoint eip;
|
|
double error = 0.0;
|
|
|
|
mesh = fes -> GetMesh();
|
|
dim = mesh -> Dimension();
|
|
e_grad.SetSize (dim);
|
|
a_grad.SetSize (dim);
|
|
Jinv.SetSize (dim);
|
|
|
|
if (norm_type & 1)
|
|
for (i = 0; i < mesh -> GetNE(); i++)
|
|
{
|
|
fe = fes -> GetFE (i);
|
|
fdof = fe -> GetDof();
|
|
transf = mesh -> GetElementTransformation (i);
|
|
el_dofs.SetSize (fdof);
|
|
dshape.SetSize (fdof, dim);
|
|
dshapet.SetSize (fdof, dim);
|
|
intorder = 2 * fe -> GetOrder(); // <----------
|
|
const IntegrationRule &ir = IntRules.Get (fe -> GetGeomType(),
|
|
intorder);
|
|
fes -> GetElementVDofs (i, vdofs);
|
|
for (k = 0; k < fdof; k++)
|
|
if (vdofs[k] >= 0)
|
|
el_dofs(k) = (*this)(vdofs[k]);
|
|
else
|
|
el_dofs(k) = - (*this)(-1-vdofs[k]);
|
|
for (j = 0; j < ir.GetNPoints(); j++)
|
|
{
|
|
const IntegrationPoint &ip = ir.IntPoint (j);
|
|
fe -> CalcDShape (ip, dshape);
|
|
transf -> SetIntPoint (&ip);
|
|
exgrad -> Eval (e_grad, *transf, ip);
|
|
CalcInverse (transf -> Jacobian(), Jinv);
|
|
Mult (dshape, Jinv, dshapet);
|
|
dshapet.MultTranspose (el_dofs, a_grad);
|
|
e_grad -= a_grad;
|
|
error += (ip.weight * transf -> Weight() *
|
|
ell_coeff -> Eval (*transf, ip) *
|
|
(e_grad * e_grad));
|
|
}
|
|
}
|
|
|
|
if (norm_type & 2)
|
|
for (i = 0; i < mesh -> GetNFaces(); i++)
|
|
{
|
|
face_elem_transf = mesh -> GetFaceElementTransformations (i);
|
|
int i1 = face_elem_transf -> Elem1No;
|
|
int i2 = face_elem_transf -> Elem2No;
|
|
intorder = fes -> GetFE (i1) -> GetOrder();
|
|
if (i2 >= 0)
|
|
if ( (k = fes -> GetFE (i2) -> GetOrder()) > intorder )
|
|
intorder = k;
|
|
intorder = 2 * intorder; // <-------------
|
|
const IntegrationRule &ir =
|
|
IntRules.Get (face_elem_transf -> FaceGeom, intorder);
|
|
err_val.SetSize (ir.GetNPoints());
|
|
ell_coeff_val.SetSize (ir.GetNPoints());
|
|
// side 1
|
|
transf = face_elem_transf -> Elem1;
|
|
fe = fes -> GetFE (i1);
|
|
fdof = fe -> GetDof();
|
|
fes -> GetElementVDofs (i1, vdofs);
|
|
shape.SetSize (fdof);
|
|
el_dofs.SetSize (fdof);
|
|
for (k = 0; k < fdof; k++)
|
|
if (vdofs[k] >= 0)
|
|
el_dofs(k) = (*this)(vdofs[k]);
|
|
else
|
|
el_dofs(k) = - (*this)(-1-vdofs[k]);
|
|
for (j = 0; j < ir.GetNPoints(); j++)
|
|
{
|
|
face_elem_transf -> Loc1.Transform (ir.IntPoint (j), eip);
|
|
fe -> CalcShape (eip, shape);
|
|
transf->SetIntPoint(&eip);
|
|
ell_coeff_val(j) = ell_coeff -> Eval (*transf, eip);
|
|
err_val(j) = exsol -> Eval (*transf, eip) - (shape * el_dofs);
|
|
}
|
|
if (i2 >= 0)
|
|
{
|
|
// side 2
|
|
transf = face_elem_transf -> Elem2;
|
|
fe = fes -> GetFE (i2);
|
|
fdof = fe -> GetDof();
|
|
fes -> GetElementVDofs (i2, vdofs);
|
|
shape.SetSize (fdof);
|
|
el_dofs.SetSize (fdof);
|
|
for (k = 0; k < fdof; k++)
|
|
if (vdofs[k] >= 0)
|
|
el_dofs(k) = (*this)(vdofs[k]);
|
|
else
|
|
el_dofs(k) = - (*this)(-1-vdofs[k]);
|
|
for (j = 0; j < ir.GetNPoints(); j++)
|
|
{
|
|
face_elem_transf -> Loc2.Transform (ir.IntPoint (j), eip);
|
|
fe -> CalcShape (eip, shape);
|
|
transf->SetIntPoint(&eip);
|
|
ell_coeff_val(j) += ell_coeff -> Eval (*transf, eip);
|
|
ell_coeff_val(j) *= 0.5;
|
|
err_val(j) -= (exsol -> Eval (*transf, eip) -
|
|
(shape * el_dofs));
|
|
}
|
|
}
|
|
transf = face_elem_transf -> Face;
|
|
for (j = 0; j < ir.GetNPoints(); j++)
|
|
{
|
|
const IntegrationPoint &ip = ir.IntPoint (j);
|
|
transf -> SetIntPoint (&ip);
|
|
error += (ip.weight * Nu * ell_coeff_val(j) *
|
|
pow (transf -> Weight(), 1.0-1.0/(dim-1)) *
|
|
err_val(j) * err_val(j));
|
|
}
|
|
}
|
|
|
|
if (error < 0.0)
|
|
return -sqrt(-error);
|
|
return sqrt(error);
|
|
}
|
|
|
|
double GridFunction::ComputeMaxError (
|
|
Coefficient *exsol[], const IntegrationRule *irs[]) const
|
|
{
|
|
double error = 0.0, a;
|
|
const FiniteElement *fe;
|
|
ElementTransformation *transf;
|
|
Vector shape;
|
|
Array<int> vdofs;
|
|
int fdof, d, i, intorder, j, k;
|
|
|
|
for (i = 0; i < fes->GetNE(); i++)
|
|
{
|
|
fe = fes->GetFE(i);
|
|
fdof = fe->GetDof();
|
|
transf = fes->GetElementTransformation(i);
|
|
shape.SetSize(fdof);
|
|
intorder = fe->GetOrder()+2; // <----------
|
|
const IntegrationRule *ir;
|
|
if (irs)
|
|
ir = irs[fe->GetGeomType()];
|
|
else
|
|
ir = &(IntRules.Get(fe->GetGeomType(), intorder));
|
|
fes->GetElementVDofs (i, vdofs);
|
|
for (j = 0; j < ir->GetNPoints(); j++)
|
|
{
|
|
const IntegrationPoint &ip = ir->IntPoint(j);
|
|
fe->CalcShape(ip, shape);
|
|
transf->SetIntPoint(&ip);
|
|
for (d = 0; d < fes->GetVDim(); d++)
|
|
{
|
|
a = 0;
|
|
for (k = 0; k < fdof; k++)
|
|
if (vdofs[fdof*d+k] >= 0)
|
|
a += (*this)(vdofs[fdof*d+k]) * shape(k);
|
|
else
|
|
a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
|
|
a -= exsol[d]->Eval(*transf, ip);
|
|
a = fabs (a);
|
|
if (error < a)
|
|
error = a;
|
|
}
|
|
}
|
|
}
|
|
|
|
return error;
|
|
}
|
|
|
|
double GridFunction::ComputeMaxError(
|
|
VectorCoefficient &exsol, const IntegrationRule *irs[]) const
|
|
{
|
|
double error = 0.0;
|
|
const FiniteElement *fe;
|
|
ElementTransformation *transf;
|
|
DenseMatrix vals, exact_vals, tr;
|
|
Vector loc_errs;
|
|
|
|
for (int i = 0; i < fes->GetNE(); i++)
|
|
{
|
|
fe = fes->GetFE(i);
|
|
transf = fes->GetElementTransformation(i);
|
|
int intorder = fe->GetOrder()+2; // <----------
|
|
const IntegrationRule *ir;
|
|
if (irs)
|
|
ir = irs[fe->GetGeomType()];
|
|
else
|
|
ir = &(IntRules.Get(fe->GetGeomType(), intorder));
|
|
GetVectorValues(i, *ir, vals, tr);
|
|
exsol.Eval(exact_vals, *transf, *ir);
|
|
vals -= exact_vals;
|
|
loc_errs.SetSize(vals.Width());
|
|
// compute the lengths of the errors at the integration points
|
|
// thus the vector max. norm is rotationally invariant
|
|
vals.Norm2(loc_errs);
|
|
double loc_error = loc_errs.Normlinf();
|
|
if (error < loc_error)
|
|
error = loc_error;
|
|
}
|
|
|
|
return error;
|
|
}
|
|
|
|
double GridFunction::ComputeW11Error (
|
|
Coefficient *exsol, VectorCoefficient *exgrad, int norm_type,
|
|
Array<int> *elems, const IntegrationRule *irs[]) const
|
|
{
|
|
// assuming vdim is 1
|
|
int i, fdof, dim, intorder, j, k;
|
|
Mesh *mesh;
|
|
const FiniteElement *fe;
|
|
ElementTransformation *transf;
|
|
Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
|
|
DenseMatrix dshape, dshapet, Jinv;
|
|
Array<int> vdofs;
|
|
double a, error = 0.0;
|
|
|
|
mesh = fes -> GetMesh();
|
|
dim = mesh -> Dimension();
|
|
e_grad.SetSize (dim);
|
|
a_grad.SetSize (dim);
|
|
Jinv.SetSize (dim);
|
|
|
|
if (norm_type & 1) // L_1 norm
|
|
for (i = 0; i < mesh -> GetNE(); i++)
|
|
{
|
|
if (elems != NULL && (*elems)[i] == 0) continue;
|
|
fe = fes -> GetFE (i);
|
|
fdof = fe -> GetDof();
|
|
transf = fes -> GetElementTransformation (i);
|
|
el_dofs.SetSize (fdof);
|
|
shape.SetSize (fdof);
|
|
intorder = fe -> GetOrder() + 1; // <----------
|
|
const IntegrationRule *ir;
|
|
if (irs)
|
|
ir = irs[fe->GetGeomType()];
|
|
else
|
|
ir = &(IntRules.Get(fe->GetGeomType(), intorder));
|
|
fes -> GetElementVDofs (i, vdofs);
|
|
for (k = 0; k < fdof; k++)
|
|
if (vdofs[k] >= 0)
|
|
el_dofs(k) = (*this)(vdofs[k]);
|
|
else
|
|
el_dofs(k) = - (*this)(-1-vdofs[k]);
|
|
for (j = 0; j < ir->GetNPoints(); j++)
|
|
{
|
|
const IntegrationPoint &ip = ir->IntPoint (j);
|
|
fe -> CalcShape (ip, shape);
|
|
transf -> SetIntPoint (&ip);
|
|
a = (el_dofs * shape) - (exsol -> Eval (*transf, ip));
|
|
error += ip.weight * transf -> Weight() * fabs (a);
|
|
}
|
|
}
|
|
|
|
if (norm_type & 2) // W^1_1 seminorm
|
|
for (i = 0; i < mesh -> GetNE(); i++)
|
|
{
|
|
if (elems != NULL && (*elems)[i] == 0) continue;
|
|
fe = fes -> GetFE (i);
|
|
fdof = fe -> GetDof();
|
|
transf = mesh -> GetElementTransformation (i);
|
|
el_dofs.SetSize (fdof);
|
|
dshape.SetSize (fdof, dim);
|
|
dshapet.SetSize (fdof, dim);
|
|
intorder = fe -> GetOrder() + 1; // <----------
|
|
const IntegrationRule *ir;
|
|
if (irs)
|
|
ir = irs[fe->GetGeomType()];
|
|
else
|
|
ir = &(IntRules.Get(fe->GetGeomType(), intorder));
|
|
fes -> GetElementVDofs (i, vdofs);
|
|
for (k = 0; k < fdof; k++)
|
|
if (vdofs[k] >= 0)
|
|
el_dofs(k) = (*this)(vdofs[k]);
|
|
else
|
|
el_dofs(k) = - (*this)(-1-vdofs[k]);
|
|
for (j = 0; j < ir->GetNPoints(); j++)
|
|
{
|
|
const IntegrationPoint &ip = ir->IntPoint (j);
|
|
fe -> CalcDShape (ip, dshape);
|
|
transf -> SetIntPoint (&ip);
|
|
exgrad -> Eval (e_grad, *transf, ip);
|
|
CalcInverse (transf -> Jacobian(), Jinv);
|
|
Mult (dshape, Jinv, dshapet);
|
|
dshapet.MultTranspose (el_dofs, a_grad);
|
|
e_grad -= a_grad;
|
|
error += ip.weight * transf -> Weight() * e_grad.Norml1();
|
|
}
|
|
}
|
|
|
|
return error;
|
|
}
|
|
|
|
double GridFunction::ComputeL1Error(
|
|
VectorCoefficient &exsol, const IntegrationRule *irs[]) const
|
|
{
|
|
double error = 0.0;
|
|
const FiniteElement *fe;
|
|
ElementTransformation *transf;
|
|
DenseMatrix vals, exact_vals, tr;
|
|
Vector loc_errs;
|
|
|
|
for (int i = 0; i < fes->GetNE(); i++)
|
|
{
|
|
fe = fes->GetFE(i);
|
|
transf = fes->GetElementTransformation(i);
|
|
int intorder = fe->GetOrder()+2; // <----------
|
|
const IntegrationRule *ir;
|
|
if (irs)
|
|
ir = irs[fe->GetGeomType()];
|
|
else
|
|
ir = &(IntRules.Get(fe->GetGeomType(), intorder));
|
|
GetVectorValues(i, *ir, vals, tr);
|
|
exsol.Eval(exact_vals, *transf, *ir);
|
|
vals -= exact_vals;
|
|
loc_errs.SetSize(vals.Width());
|
|
// compute the lengths of the errors at the integration points
|
|
// thus the vector L_1 norm is rotationally invariant
|
|
vals.Norm2(loc_errs);
|
|
for (int j = 0; j < ir->GetNPoints(); j++)
|
|
{
|
|
const IntegrationPoint &ip = ir->IntPoint (j);
|
|
transf->SetIntPoint(&ip);
|
|
error += ip.weight * transf->Weight() * loc_errs(j);
|
|
}
|
|
}
|
|
|
|
return error;
|
|
}
|
|
|
|
GridFunction & GridFunction::operator= (double value)
|
|
{
|
|
for(int i = 0; i < size; i++)
|
|
data[i] = value;
|
|
return *this;
|
|
}
|
|
|
|
GridFunction & GridFunction::operator= (const Vector &v)
|
|
{
|
|
for(int i = 0; i < size; i++)
|
|
data[i] = v(i);
|
|
return *this;
|
|
}
|
|
|
|
GridFunction & GridFunction::operator= (const GridFunction &v)
|
|
{
|
|
return this->operator=((const Vector &)v);
|
|
}
|
|
|
|
void GridFunction::Save (ostream &out)
|
|
{
|
|
fes -> Save (out);
|
|
out << endl;
|
|
Vector::Print (out, 1);
|
|
}
|
|
|
|
void GridFunction::SaveVTK(ostream &out, const string &field_name, int ref)
|
|
{
|
|
Mesh *mesh = fes->GetMesh();
|
|
RefinedGeometry *RefG;
|
|
Vector val;
|
|
DenseMatrix vval, pmat;
|
|
|
|
if (VectorDim() == 1)
|
|
{
|
|
// scalar data
|
|
out << "SCALARS " << field_name << " double 1\n"
|
|
<< "LOOKUP_TABLE default\n";
|
|
for (int i = 0; i < mesh->GetNE(); i++)
|
|
{
|
|
RefG = GlobGeometryRefiner.Refine(
|
|
mesh->GetElementBaseGeometry(i), ref, 1);
|
|
|
|
GetValues(i, RefG->RefPts, val, pmat);
|
|
|
|
for (int j = 0; j < val.Size(); j++)
|
|
{
|
|
out << val(j) << '\n';
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// vector data
|
|
out << "VECTORS " << field_name << " double\n";
|
|
for (int i = 0; i < mesh->GetNE(); i++)
|
|
{
|
|
RefG = GlobGeometryRefiner.Refine(
|
|
mesh->GetElementBaseGeometry(i), ref, 1);
|
|
|
|
GetVectorValues(i, RefG->RefPts, vval, pmat);
|
|
|
|
for (int j = 0; j < vval.Width(); j++)
|
|
{
|
|
out << vval(0, j) << ' ' << vval(1, j) << ' ';
|
|
if (vval.Height() == 2)
|
|
out << 0.0;
|
|
else
|
|
out << vval(2, j);
|
|
out << '\n';
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void GridFunction::SaveSTLTri (ostream &out, double p1[], double p2[],
|
|
double p3[])
|
|
{
|
|
double v1[3] = { p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2] };
|
|
double v2[3] = { p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2] };
|
|
double n[] = { v1[1] * v2[2] - v1[2] * v2[1],
|
|
v1[2] * v2[0] - v1[0] * v2[2],
|
|
v1[0] * v2[1] - v1[1] * v2[0] };
|
|
double rl = 1.0 / sqrt (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
|
|
n[0] *= rl; n[1] *= rl; n[2] *= rl;
|
|
|
|
out << " facet normal " << n[0] << ' ' << n[1] << ' ' << n[2]
|
|
<< "\n outer loop"
|
|
<< "\n vertex " << p1[0] << ' ' << p1[1] << ' ' << p1[2]
|
|
<< "\n vertex " << p2[0] << ' ' << p2[1] << ' ' << p2[2]
|
|
<< "\n vertex " << p3[0] << ' ' << p3[1] << ' ' << p3[2]
|
|
<< "\n endloop\n endfacet\n";
|
|
}
|
|
|
|
void GridFunction::SaveSTL (ostream &out, int TimesToRefine)
|
|
{
|
|
Mesh *mesh = fes -> GetMesh();
|
|
|
|
if (mesh -> Dimension() != 2)
|
|
return;
|
|
|
|
int i, j, k, l, n;
|
|
DenseMatrix pointmat;
|
|
Vector values;
|
|
RefinedGeometry * RefG;
|
|
double pts[4][3], bbox[3][2];
|
|
|
|
out << "solid GridFunction\n";
|
|
|
|
bbox[0][0] = bbox[0][1] = bbox[1][0] = bbox[1][1] =
|
|
bbox[2][0] = bbox[2][1] = 0.0;
|
|
for (i = 0; i < mesh -> GetNE(); i++)
|
|
{
|
|
n = fes -> GetFE(i) -> GetGeomType();
|
|
RefG = GlobGeometryRefiner.Refine (n, TimesToRefine);
|
|
GetValues (i, RefG -> RefPts, values, pointmat);
|
|
Array<int> &RG = RefG -> RefGeoms;
|
|
n = Geometries.NumBdr (n);
|
|
for (k = 0; k < RG.Size()/n; k++)
|
|
{
|
|
for (j = 0; j < n; j++)
|
|
{
|
|
l = RG[n*k+j];
|
|
pts[j][0] = pointmat(0,l);
|
|
pts[j][1] = pointmat(1,l);
|
|
pts[j][2] = values(l);
|
|
}
|
|
|
|
if (n == 3)
|
|
{
|
|
SaveSTLTri (out, pts[0], pts[1], pts[2]);
|
|
}
|
|
else
|
|
{
|
|
SaveSTLTri (out, pts[0], pts[1], pts[2]);
|
|
SaveSTLTri (out, pts[0], pts[2], pts[3]);
|
|
}
|
|
}
|
|
|
|
if (i == 0)
|
|
{
|
|
bbox[0][0] = pointmat(0,0);
|
|
bbox[0][1] = pointmat(0,0);
|
|
bbox[1][0] = pointmat(1,0);
|
|
bbox[1][1] = pointmat(1,0);
|
|
bbox[2][0] = values(0);
|
|
bbox[2][1] = values(0);
|
|
}
|
|
|
|
for (j = 0; j < values.Size(); j++)
|
|
{
|
|
if (bbox[0][0] > pointmat(0,j))
|
|
bbox[0][0] = pointmat(0,j);
|
|
if (bbox[0][1] < pointmat(0,j))
|
|
bbox[0][1] = pointmat(0,j);
|
|
if (bbox[1][0] > pointmat(1,j))
|
|
bbox[1][0] = pointmat(1,j);
|
|
if (bbox[1][1] < pointmat(1,j))
|
|
bbox[1][1] = pointmat(1,j);
|
|
if (bbox[2][0] > values(j))
|
|
bbox[2][0] = values(j);
|
|
if (bbox[2][1] < values(j))
|
|
bbox[2][1] = values(j);
|
|
}
|
|
}
|
|
|
|
cout << "[xmin,xmax] = [" << bbox[0][0] << ',' << bbox[0][1] << "]\n"
|
|
<< "[ymin,ymax] = [" << bbox[1][0] << ',' << bbox[1][1] << "]\n"
|
|
<< "[zmin,zmax] = [" << bbox[2][0] << ',' << bbox[2][1] << ']'
|
|
<< endl;
|
|
|
|
out << "endsolid GridFunction" << endl;
|
|
}
|
|
|
|
|
|
void ComputeFlux(BilinearFormIntegrator &blfi,
|
|
GridFunction &u,
|
|
GridFunction &flux, int wcoef, int sd)
|
|
{
|
|
int i, j, nfe;
|
|
FiniteElementSpace *ufes, *ffes;
|
|
ElementTransformation *Transf;
|
|
|
|
ufes = u.FESpace();
|
|
ffes = flux.FESpace();
|
|
nfe = ufes->GetNE();
|
|
Array<int> udofs;
|
|
Array<int> fdofs;
|
|
Array<int> overlap(flux.Size());
|
|
Vector ul, fl;
|
|
|
|
flux = 0.0;
|
|
|
|
for (i = 0; i < overlap.Size(); i++)
|
|
overlap[i] = 0;
|
|
|
|
for (i = 0; i < nfe; i++)
|
|
if (sd < 0 || ufes->GetAttribute(i) == sd)
|
|
{
|
|
ufes->GetElementVDofs(i, udofs);
|
|
ffes->GetElementVDofs(i, fdofs);
|
|
|
|
ul.SetSize(udofs.Size());
|
|
for (j = 0; j < ul.Size(); j++)
|
|
ul(j) = u(udofs[j]);
|
|
|
|
Transf = ufes->GetElementTransformation(i);
|
|
blfi.ComputeElementFlux(*ufes->GetFE(i), *Transf, ul,
|
|
*ffes->GetFE(i), fl, wcoef);
|
|
|
|
flux.AddElementVector(fdofs, fl);
|
|
|
|
for (j = 0; j < fdofs.Size(); j++)
|
|
overlap[fdofs[j]]++;
|
|
}
|
|
|
|
for (i = 0; i < overlap.Size(); i++)
|
|
if (overlap[i] != 0)
|
|
flux(i) /= overlap[i];
|
|
}
|
|
|
|
void ZZErrorEstimator(BilinearFormIntegrator &blfi,
|
|
GridFunction &u,
|
|
GridFunction &flux, Vector &ErrorEstimates,
|
|
int wsd)
|
|
{
|
|
int i, j, s, nfe, nsd;
|
|
FiniteElementSpace *ufes, *ffes;
|
|
ElementTransformation *Transf;
|
|
|
|
ufes = u.FESpace();
|
|
ffes = flux.FESpace();
|
|
nfe = ufes->GetNE();
|
|
Array<int> udofs;
|
|
Array<int> fdofs;
|
|
Vector ul, fl, fla;
|
|
|
|
ErrorEstimates.SetSize(nfe);
|
|
|
|
nsd = 1;
|
|
if (wsd)
|
|
for (i = 0; i < nfe; i++)
|
|
if ( (j=ufes->GetAttribute(i)) > nsd)
|
|
nsd = j;
|
|
|
|
for (s = 1; s <= nsd; s++)
|
|
{
|
|
if (wsd)
|
|
ComputeFlux(blfi, u, flux, 0, s);
|
|
else
|
|
ComputeFlux(blfi, u, flux, 0);
|
|
|
|
for (i = 0; i < nfe; i++)
|
|
if (!wsd || ufes->GetAttribute(i) == s)
|
|
{
|
|
ufes->GetElementVDofs(i, udofs);
|
|
ffes->GetElementVDofs(i, fdofs);
|
|
|
|
ul.SetSize(udofs.Size());
|
|
for (j = 0; j < ul.Size(); j++)
|
|
ul(j) = u(udofs[j]);
|
|
|
|
fla.SetSize(fdofs.Size());
|
|
for (j = 0; j < fla.Size(); j++)
|
|
fla(j) = flux(fdofs[j]);
|
|
|
|
Transf = ufes->GetElementTransformation(i);
|
|
blfi.ComputeElementFlux(*ufes->GetFE(i), *Transf, ul,
|
|
*ffes->GetFE(i), fl, 0);
|
|
|
|
fl -= fla;
|
|
|
|
ErrorEstimates(i) = blfi.ComputeFluxEnergy(*ffes->GetFE(i),
|
|
*Transf, fl);
|
|
}
|
|
}
|
|
}
|