Permalink
Browse files

Merge branch 'quadsurface' into 'master'

Quadsurface

See merge request jschoeberl/ngsolve!293
  • Loading branch information...
JSchoeberl committed Jan 30, 2018
2 parents 81d85c9 + 1b8c272 commit a3771b90c544a102e48749ee83eaeefa27668762
Showing with 52 additions and 20 deletions.
  1. +19 −4 comp/facetsurffespace.cpp
  2. +3 −0 comp/facetsurffespace.hpp
  3. +14 −9 comp/hdivdivsurfacespace.cpp
  4. +5 −6 comp/hdivhosurfacefespace.cpp
  5. +11 −1 python/utils.py
@@ -243,6 +243,18 @@ namespace ngcomp
}
}
template <ELEMENT_TYPE ET>
FiniteElement & FacetSurfaceFESpace :: T_GetFE (int elnr, Allocator & alloc) const
{
Ngs_Element ngel = ma->GetElement<ET_trait<ET>::DIM,VOL> (elnr);
FacetFE<ET> * fe = new (alloc) FacetFE<ET> ();
fe -> SetVertexNumbers (ngel.Vertices());
fe -> SetOrder (order);
fe -> ComputeNDof();
return *fe;
}
// ------------------------------------------------------------------------
FiniteElement & FacetSurfaceFESpace :: GetFE (ElementId ei, Allocator & lh) const
@@ -258,19 +270,22 @@ namespace ngcomp
}
case BND:
{
FacetFE<ET_TRIG>* fe = 0;
//FacetFE<ET_TRIG>* fet = 0;
//FacetFE<ET_QUAD>* feq = 0;
switch (ma->GetElType(ei))
{
case ET_TRIG: fe = new (lh) FacetFE<ET_TRIG> (); break;
case ET_TRIG: return T_GetFE<ET_TRIG>(ei.Nr(), lh);//fe = new (lh) FacetFE<ET_TRIG> (); break;
case ET_QUAD: return T_GetFE<ET_QUAD>(ei.Nr(), lh);//fe = new (lh) FacetFE<ET_QUAD> (); break;
default:
throw Exception (string("FacetSurfaceFESpace::GetFE: unsupported element ")+
ElementTopology::GetElementName(ma->GetElType(ei)));
}
switch (ma->GetElType(ei))
/*switch (ma->GetElType(ei))
{
case ET_TRIG:
case ET_QUAD:
{
fe -> SetVertexNumbers (vnums);
fe -> SetOrder (order);
@@ -281,7 +296,7 @@ namespace ngcomp
default:
break;
}
return *fe;
return *fe;*/
}
case BBND:
{
@@ -46,6 +46,9 @@ namespace ngcomp
///
virtual size_t GetNDofLevel (int level) const override;
template <ELEMENT_TYPE ET>
FiniteElement & T_GetFE (int elnr, Allocator & alloc) const;
virtual FiniteElement & GetFE (ElementId ei, Allocator & lh) const override;
using FESpace::GetDofNrs;
@@ -281,7 +281,7 @@ namespace ngcomp
{
first_face_dof[i] = ndof;
ma->GetEdgePNums(i, pnums);
pnums = ma->GetEdgePNums(i);
switch (pnums.Size())
{
@@ -327,7 +327,7 @@ namespace ngcomp
if (noncontinuous)
{
cout << "Update before discont" << endl;
//cout << "Update before discont" << endl;
ndof = 0;
Array<int> pnums;
@@ -359,8 +359,8 @@ namespace ngcomp
}
first_element_dof[nel] = ndof;
first_face_dof = 0;
cout << "fed" << first_element_dof << endl;
cout << "ffd" << first_face_dof << endl;
//cout << "fed" << first_element_dof << endl;
//cout << "ffd" << first_face_dof << endl;
}
@@ -417,11 +417,16 @@ namespace ngcomp
fe = trigfe;
break;
}
// case ET_QUAD:
// fe = new (lh.Alloc(sizeof(HDivSymQuad))) HDivSymQuad();
// break;
default:
cerr << "element type " << int(ma->GetElType(ei)) << " not there in hdivsymsurf" << endl;
case ET_QUAD:
{
auto * quadfe = new (lh) HDivDivFE<ET_QUAD> (order, false /* plus */);
quadfe->SetVertexNumbers(vnums);
quadfe->ComputeNDof();
fe = quadfe;
break;
}
default:
cerr << "element type " << int(ma->GetElType(ei)) << " not there in hdivdivsurf" << endl;
}
ArrayMem<INT<2>,4> order_ed;
@@ -456,12 +456,11 @@ class DiffOpDivHDivSurface : public DiffOp<DiffOpDivHDivSurface<D, FEL> >
//inci = pc[0]*(pc[0]-1)/2;
break;
case ET_QUAD:
throw Exception("not implemented yet");
//if (!ho_div_free)
// inci = pc[0]*pc[1] + p[0]*p[1]+p[0]+p[1];
//else
// inci = pc[0]*pc[1];
//break;
if (!ho_div_free)
inci = p[0]*p[1] + p[0]*p[1]+p[0]+p[1];
else
inci = p[0]*p[1];
break;
default: // for the compiler
break;
}
@@ -131,8 +131,18 @@ def Cof(m):
-m[0,0]*m[1,2]+m[1,0]*m[0,2],
m[0,0]*m[1,1]-m[1,0]*m[0,1] ), dims=(3,3) )
def Inv(m):
return 1/Det(m)*Cof(m).trans
def Sym(m):
return 0.5*(m+m.trans)
__all__ = ['x', 'y', 'z', 'Laplace', 'Mass', 'Source', 'Neumann', 'H1', 'VectorH1', 'FacetFESpace', 'L2', 'VectorL2', 'SurfaceL2', 'HDivDiv', 'NumberSpace', 'grad', 'curl', 'div','Mesh', 'ConstantCF', 'DomainConstantCF', 'Id', 'Trace', 'Det', 'Cross', 'Cof']
def Skew(m):
return 0.5*(m-m.trans)
def OuterProduct(a, b):
return CoefficientFunction( tuple([a[i]*b[j] for i in range(a.dim) for j in range(b.dim)]), dims=(a.dim,b.dim) )
__all__ = ['x', 'y', 'z', 'Laplace', 'Mass', 'Source', 'Neumann', 'H1', 'VectorH1', 'FacetFESpace', 'L2', 'VectorL2', 'SurfaceL2', 'HDivDiv', 'NumberSpace', 'grad', 'curl', 'div','Mesh', 'ConstantCF', 'DomainConstantCF', 'Id', 'Trace', 'Det', 'Cross', 'Cof', 'Inv', 'Sym', 'Skew', 'OuterProduct']

0 comments on commit a3771b9

Please sign in to comment.