From e86a825bee472f48118667abb95057cacd951d2a Mon Sep 17 00:00:00 2001
From: Boris Sedji <sedji.boris@hotmail.com>
Date: Wed, 3 Feb 2010 15:55:24 +0000
Subject: [PATCH]
---
contrib/arc/Classes/FuncGradDisc.h | 192 +++++++++
contrib/arc/Classes/OctreeSolver.cpp | 488 ++++++++++++++++++++++
contrib/arc/Classes/OctreeSolver.h | 53 +++
contrib/arc/Classes/gLevelSetMesh.cpp | 21 +
contrib/arc/Classes/gLevelSetMesh.h | 61 +++
contrib/arc/Classes/xFemFunctionSpace.cpp | 220 ++++++++++
contrib/arc/Classes/xFemFunctionSpace.h | 111 +++++
7 files changed, 1146 insertions(+)
create mode 100644 contrib/arc/Classes/FuncGradDisc.h
create mode 100644 contrib/arc/Classes/OctreeSolver.cpp
create mode 100644 contrib/arc/Classes/OctreeSolver.h
create mode 100644 contrib/arc/Classes/gLevelSetMesh.cpp
create mode 100644 contrib/arc/Classes/gLevelSetMesh.h
create mode 100644 contrib/arc/Classes/xFemFunctionSpace.cpp
create mode 100644 contrib/arc/Classes/xFemFunctionSpace.h
diff --git a/contrib/arc/Classes/FuncGradDisc.h b/contrib/arc/Classes/FuncGradDisc.h
new file mode 100644
index 0000000000..7bfcde4631
--- /dev/null
+++ b/contrib/arc/Classes/FuncGradDisc.h
@@ -0,0 +1,192 @@
+//
+// Description : Heaviside function based on level set discontinuity description
+//
+//
+// Author: <Boris Sedji>, 01/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef _FUNCGRADDISC_H_
+#define _FUNCGRADDISC_H_
+
+#include "simpleFunction.h"
+#include "DILevelset.h"
+#include "MVertex.h"
+#include "GModel.h"
+//#include "gLevelSetMesh.cpp"
+
+
+class FuncGradDisc : public simpleFunction<double> {
+ private :
+
+ gLevelset *_ls;
+ GModel * _pModel;
+
+
+ public :
+ FuncGradDisc(gLevelset *ls,GModel *pModel){
+ _ls = ls;
+ _pModel = pModel;
+
+ }
+ virtual double operator () (double x,double y,double z) const
+ {
+ double test = (*_ls)(x,y,z);
+ SPoint3 p(x,y,z);
+ MElement *e = _pModel->getMeshElementByCoord(p);
+ if (e->getParent()) e = e->getParent();
+ int numV = e->getNumVertices();
+ double f = 0;
+ for (int i = 0;i<numV;i++)
+ {
+ f = f + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()));
+ }
+ f = f - abs((*_ls)(x,y,z));
+ return f;
+ }
+
+ double operator () (double x,double y,double z, MElement *e) const {
+
+ SPoint3 p(x,y,z);
+ if (e->getParent()) e = e->getParent();
+ int numV = e->getNumVertices();
+ double xyz[3] = {x,y,z};
+ double uvw[3];
+ e->xyz2uvw(xyz,uvw);
+ double val[30];
+ e->getShapeFunctions(uvw[0], uvw[1], uvw[2], val);
+ double f = 0;
+ for (int i = 0;i<numV;i++)
+ {
+ f = f + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*val[i];
+ }
+ f = f - abs((*_ls)(x,y,z));
+ return f;
+
+ }
+
+ virtual void gradient (double x, double y, double z,
+ double & dfdx, double & dfdy, double & dfdz) const
+ {
+
+ SPoint3 p(x,y,z);
+ MElement *e = _pModel->getMeshElementByCoord(p);
+ if (e->getParent()) e = e->getParent();
+ int numV = e->getNumVertices();
+ double xyz[3] = {x,y,z};
+ double uvw[3];
+ e->xyz2uvw(xyz,uvw);
+ double gradsuvw[256][3];
+ e->getGradShapeFunctions(uvw[0],uvw[1],uvw[2],gradsuvw);
+
+ double jac[3][3];
+ double invjac[3][3];
+ double dNdx;
+ double dNdy;
+ double dNdz;
+ const double detJ = e->getJacobian(uvw[0], uvw[1], uvw[2], jac);
+ inv3x3(jac, invjac);
+
+ dfdx = 0;
+ dfdy = 0;
+ dfdz = 0;
+
+ if ((*_ls)(x,y,z)>0)
+ {
+ for (int i = 0;i<numV;i++)
+ {
+ dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
+ dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
+ dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
+
+ dfdx = dfdx + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdx ;
+ dfdx = dfdx - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdx;
+ dfdy = dfdy + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdy ;
+ dfdy = dfdy - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdy;
+ dfdz = dfdz + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdz ;
+ dfdz = dfdz - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdz;
+ }
+ }else
+ {
+ for (int i = 0;i<numV;i++)
+ {
+ dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
+ dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
+ dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
+
+ dfdx = dfdx + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdx ;
+ dfdx = dfdx + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdx;
+ dfdy = dfdy + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdy ;
+ dfdy = dfdy + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdy;
+ dfdz = dfdz + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdz ;
+ dfdz = dfdz + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdz;
+ }
+ }
+
+
+
+ }
+
+
+ void gradient (double x, double y, double z,
+ double & dfdx, double & dfdy, double & dfdz,MElement *e) const
+ {
+
+ SPoint3 p(x,y,z);
+ if (e->getParent()) e = e->getParent();
+ int numV = e->getNumVertices();
+ double xyz[3] = {x,y,z};
+ double uvw[3];
+ e->xyz2uvw(xyz,uvw);
+ double gradsuvw[256][3];
+ e->getGradShapeFunctions(uvw[0],uvw[1],uvw[2],gradsuvw);
+
+ double jac[3][3];
+ double invjac[3][3];
+ double dNdx;
+ double dNdy;
+ double dNdz;
+ const double detJ = e->getJacobian(uvw[0], uvw[1], uvw[2], jac);
+ inv3x3(jac, invjac);
+
+ dfdx = 0;
+ dfdy = 0;
+ dfdz = 0;
+
+ if ((*_ls)(x,y,z)<0)
+ {
+ for (int i = 0;i<numV;i++)
+ {
+ dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
+ dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
+ dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
+
+ dfdx = dfdx + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdx ;
+ dfdx = dfdx - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdx;
+ dfdy = dfdy + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdy ;
+ dfdy = dfdy - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdy;
+ dfdz = dfdz + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdz ;
+ dfdz = dfdz - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdz;
+ }
+ }else
+ {
+ for (int i = 0;i<numV;i++)
+ {
+ dNdx = invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2];
+ dNdy = invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2];
+ dNdz = invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2];
+
+ dfdx = dfdx + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdx ;
+ dfdx = dfdx + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdx;
+ dfdy = dfdy + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdy ;
+ dfdy = dfdy + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdy;
+ dfdz = dfdz + abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*dNdz ;
+ dfdz = dfdz + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdz;
+ }
+ }
+ }
+};
+
+#endif
diff --git a/contrib/arc/Classes/OctreeSolver.cpp b/contrib/arc/Classes/OctreeSolver.cpp
new file mode 100644
index 0000000000..2401e6f013
--- /dev/null
+++ b/contrib/arc/Classes/OctreeSolver.cpp
@@ -0,0 +1,488 @@
+//
+// Description :
+//
+//
+// Author: <Boris Sedji>, 01/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+
+
+#include <string.h>
+#include "GmshConfig.h"
+#include "elasticitySolver.h"
+#include "linearSystemCSR.h"
+#include "linearSystemPETSc.h"
+#include "linearSystemGMM.h"
+#include "Numeric.h"
+#include "GModel.h"
+#include "functionSpace.h"
+#include "terms.h"
+#include "solverAlgorithms.h"
+#include "quadratureRules.h"
+#include "solverField.h"
+#if defined(HAVE_POST)
+#include "PView.h"
+#include "PViewData.h"
+#include "MPoint.h"
+#endif
+
+#include "OctreeSolver.h"
+#include "DILevelset.h"
+#include "NodeEnrichedFS.cpp"
+#include "gLevelSetMesh.h"
+
+#include "groupOfElements.h"
+#include "FuncGradDisc.h"
+
+#include "xFemFunctionSpace.cpp"
+
+
+OctreeSolver::~OctreeSolver()
+{
+ //delete _funcEnrichment;
+}
+
+
+void OctreeSolver::setMesh(const std::string &meshFileName)
+{
+ std::cout<< "Mesh file : " << meshFileName.c_str() <<"\n";
+ GModel *p = new GModel();
+ p->readMSH(meshFileName.c_str());
+ _dim = p->getNumRegions() ? 3 : 2;
+ int tag = 1;
+ //if (_ls) delete _ls;
+ gLevelSetMesh *ls = new gLevelSetMesh(&_LevelSetValue,p,tag);
+ GModel *pcut = p->buildCutGModel(ls);
+ pModel = pcut;
+ std::string ModelName = meshFileName.c_str() ;
+ ModelName.erase(ModelName.find(".",0),4);
+ ModelName.insert(ModelName.length(),"_cutLS.msh");
+ pcut->writeMSH(ModelName,2.1,false,false);
+ GModel::setCurrent(pModel);
+ _ls = new gLevelSetMesh(&_LevelSetValue,p,tag);
+ //delete p;
+ delete ls;
+}
+
+void OctreeSolver::setLevelSetValue(const std::string &fn)
+{
+ FILE *f = fopen(fn.c_str(), "r");
+ char what[256];
+ int num;
+ float val;
+ int numnodes;
+ fscanf(f, "%s", what);
+ std::cout<<" Fichier levelset :"<< fn << "\n";
+ printf("%s",what);
+ fscanf(f, "%d", &numnodes);
+ std::cout<<" numnodes :"<< numnodes << "\n";
+ while(!feof(f)){
+ fscanf(f,"%d %f",&num,&val);
+ _LevelSetValue[num] = val ;
+ //std::cout<<"_LevelSetValue : "<< val << "\n";
+ }
+ fclose(f);
+}
+
+
+
+
+void OctreeSolver::readInputFile(const std::string &fn)
+{
+ FILE *f = fopen(fn.c_str(), "r");
+ char what[256];
+ while(!feof(f)){
+ if(fscanf(f, "%s", what) != 1) return;
+ if (!strcmp(what, "ElasticDomain")){
+ GModel::setCurrent(pModel);
+ elasticField field;
+ int physical;
+ if(fscanf(f, "%d %lf %lf", &physical, &field._E, &field._nu) != 3) return;
+ field._tag = _tag;
+ field.g = new groupOfElements (_dim, physical);
+ elasticFields.push_back(field);
+ }
+ else if (!strcmp(what, "Void")){
+ // elasticField field;
+ // create enrichment ...
+ // create the group ...
+ // assign a tag
+ // elasticFields.push_back(field);
+ }
+ else if (!strcmp(what, "NodeDisplacement")){
+ double val;
+ int node, comp;
+ if(fscanf(f, "%d %d %lf", &node, &comp, &val) != 3) return;
+ dirichletBC diri;
+ diri.g = new groupOfElements (0, node);
+ diri._f= simpleFunction<double>(val);
+ diri._comp=comp;
+ diri._tag=node;
+ diri.onWhat=BoundaryCondition::ON_VERTEX;
+ allDirichlet.push_back(diri);
+ }
+ else if (!strcmp(what, "EdgeDisplacement")){
+ GModel::setCurrent(pModel);
+ double val;
+ int edge, comp;
+ if(fscanf(f, "%d %d %lf", &edge, &comp, &val) != 3) return;
+ dirichletBC diri;
+ diri.g = new groupOfElements (1, edge);
+ diri._f= simpleFunction<double>(val);
+ diri._comp=comp;
+ diri._tag=edge;
+ diri.onWhat=BoundaryCondition::ON_EDGE;
+ allDirichlet.push_back(diri);
+ }
+ else if (!strcmp(what, "FaceDisplacement")){
+ double val;
+ int face, comp;
+ if(fscanf(f, "%d %d %lf", &face, &comp, &val) != 3) return;
+ dirichletBC diri;
+ diri.g = new groupOfElements (2, face);
+ diri._f= simpleFunction<double>(val);
+ diri._comp=comp;
+ diri._tag=face;
+ diri.onWhat=BoundaryCondition::ON_FACE;
+ allDirichlet.push_back(diri);
+ }
+ else if (!strcmp(what, "NodeForce")){
+ double val1, val2, val3;
+ int node;
+ if(fscanf(f, "%d %lf %lf %lf", &node, &val1, &val2, &val3) != 4) return;
+ neumannBC neu;
+ neu.g = new groupOfElements (0, node);
+ neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+ neu._tag=node;
+ neu.onWhat=BoundaryCondition::ON_VERTEX;
+ allNeumann.push_back(neu);
+ }
+ else if (!strcmp(what, "EdgeForce")){
+ GModel::setCurrent(pModel);
+ double val1, val2, val3;
+ int edge;
+ if(fscanf(f, "%d %lf %lf %lf", &edge, &val1, &val2, &val3) != 4) return;
+ neumannBC neu;
+ neu.g = new groupOfElements (1, edge);
+ neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+ neu._tag=edge;
+ neu.onWhat=BoundaryCondition::ON_EDGE;
+ allNeumann.push_back(neu);
+ }
+ else if (!strcmp(what, "FaceForce")){
+ double val1, val2, val3;
+ int face;
+ if(fscanf(f, "%d %lf %lf %lf", &face, &val1, &val2, &val3) != 4) return;
+ neumannBC neu;
+ neu.g = new groupOfElements (2, face);
+ neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+ neu._tag=face;
+ neu.onWhat=BoundaryCondition::ON_FACE;
+ allNeumann.push_back(neu);
+ }
+ else if (!strcmp(what, "VolumeForce")){
+ double val1, val2, val3;
+ int volume;
+ if(fscanf(f, "%d %lf %lf %lf", &volume, &val1, &val2, &val3) != 4) return;
+ neumannBC neu;
+ neu.g = new groupOfElements (3, volume);
+ neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+ neu._tag=volume;
+ neu.onWhat=BoundaryCondition::ON_VOLUME;
+ allNeumann.push_back(neu);
+ }
+ else if (!strcmp(what, "LevelSetFile")){
+ char name[245];
+ if(fscanf(f, "%s", name) != 1) return;
+ setLevelSetValue(name);
+ }
+ else if (!strcmp(what, "MeshFile")){
+ char name[245];
+ if(fscanf(f, "%s", name) != 1) return;
+ setMesh(name);
+ }
+ else {
+ Msg::Error("Invalid input : %s", what);
+// return;
+ }
+ }
+ fclose(f);
+}
+
+
+void OctreeSolver::solve(){
+
+ //GModel::setCurrent(pModel);
+
+#if defined(HAVE_TAUCS)
+ linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+#elif defined(HAVE_PETSC)
+ linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
+#else
+ linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
+ lsys->setNoisy(2);
+#endif
+ if (pAssembler) delete pAssembler;
+ pAssembler = new dofManager<double>(lsys);
+
+ // determine all enriched nodes and save in map member _TagEnrichedVertex
+
+ for (int i = 0; i < elasticFields.size(); ++i){
+ std::set<MElement*>::const_iterator it = elasticFields[i].g->begin();
+ for ( ; it != elasticFields[i].g->end(); ++it)
+ {
+ MElement *e = *it;
+ if (e->getParent()) // if element got parents
+ {
+ for (int k = 0; k < e->getParent()->getNumVertices(); ++k)
+ { // for all vertices in the element parent
+ _TagEnrichedVertex.insert(e->getParent()->getVertex(k)->getNum());
+ }
+ }
+ }
+ }
+
+ //_TagEnrichedVertex.clear();
+ _EnrichComp.insert(0);
+ _EnrichComp.insert(1);
+ _EnrichComp.insert(2);
+
+ _funcEnrichment = new FuncGradDisc(_ls,pModel);
+
+
+ // space function definition //
+ // Lagrange space with normal dofs
+ FunctionSpace<SVector3> *LagrangeSpace;
+ if (_dim==3) LagrangeSpace=new VectorLagrangeFunctionSpace(_tag);
+ if (_dim==2) LagrangeSpace=new VectorLagrangeFunctionSpace(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+
+ // Filtered space
+ FilterNodeEnriched *filter = new FilterNodeEnriched(&_TagEnrichedVertex , &_EnrichComp);
+ FilteredFS<SVector3> *LagrangeFiltered = new FilteredFS<SVector3>(LagrangeSpace,filter);
+
+ // xfem filtered space
+ xFemFS<SVector3> *xFemLagrange = new xFemFS<SVector3>(LagrangeFiltered,_funcEnrichment);
+
+ // Composite space : xfem + Lagrange
+ if (LagSpace) delete LagSpace;
+
+ LagSpace = new CompositeFunctionSpace<SVector3>(*LagrangeSpace, *xFemLagrange);
+
+
+ std::cout << "Dirichlet BC"<< std::endl;
+ for (unsigned int i = 0; i < allDirichlet.size(); i++)
+ {
+ FilterDofComponent filter(allDirichlet[i]._comp);
+ FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,allDirichlet[i]._f,filter);
+ }
+
+ // we number the dofs : when a dof is numbered, it cannot be numbered
+ // again with another number.
+ for (unsigned int i = 0; i < elasticFields.size(); ++i)
+ {
+ NumberDofs(*LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),*pAssembler);
+ }
+
+ // Now we start the assembly process
+ // First build the force vector
+
+ GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
+ std::cout << "Neumann BC"<< std::endl;
+ for (unsigned int i = 0; i < allNeumann.size(); i++)
+ {
+ LoadTerm<SVector3> Lterm(*LagSpace,allNeumann[i]._f);
+ Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler);
+ }
+
+// bulk material law
+
+ GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
+ for (unsigned int i = 0; i < elasticFields.size(); i++)
+ {
+ IsotropicElasticTerm Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu);
+// LaplaceTerm<SVector3,SVector3> Eterm(*LagSpace);
+ Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
+ }
+
+ printf("-- done assembling!\n");
+ std::cout<<"Dof Number : " << pAssembler->sizeOfR() <<"\n";
+ lsys->systemSolve();
+ printf("-- done solving!\n");
+ double energ=0;
+ for (unsigned int i = 0; i < elasticFields.size(); i++)
+ {
+ SolverField<SVector3> Field(pAssembler, LagSpace);
+ IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+ BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+ Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
+ }
+ printf("elastic energy=%f\n",energ);
+
+//
+// // determine all enriched nodes and save in map member _TagEnrichedVertex
+//
+//// for (int i = 0; i < elasticFields.size(); ++i){
+//// std::set<MElement*>::const_iterator it = elasticFields[i].g->begin();
+//// for ( ; it != elasticFields[i].g->end(); ++it)
+//// {
+//// MElement *e = *it;
+//// if (e->getParent()) // if element got parents
+//// {
+//// for (int k = 0; k < e->getParent()->getNumVertices(); ++k)
+//// { // for all vertices in the element parent
+//// _TagEnrichedVertex.insert(e->getParent()->getVertex(k)->getNum());
+//// }
+//// }
+//// }
+//// }
+//
+// _TagEnrichedVertex.clear();
+//
+// // enriched composant
+// _EnrichComp.push_back(0);
+// _EnrichComp.push_back(1);
+// _EnrichComp.push_back(2);
+//
+// // level set definition (in .dat ??)
+// double a(0.), b(1.), c(0.), d(-4.7);
+// int n(1); // pour level set
+// gLevelset *ls = new gLevelsetPlane(a, b, c, d, n);
+// _funcEnrichment = new FuncHeaviside(ls);
+//
+//
+// // space function definition
+// if (LagSpace) delete LagSpace;
+// FunctionSpace<SVector3> *NLagSpace;
+// if (LagSpace) delete LagSpace;
+// if (_dim==3) NLagSpace=new VectorLagrangeFunctionSpace(_tag);
+// if (_dim==2) NLagSpace=new VectorLagrangeFunctionSpace(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+// LagSpace = NLagSpace;
+// //LagSpace = new NodeEnrichedFS<SVector3>(NLagSpace, &_TagEnrichedVertex ,&_EnrichComp, _funcEnrichment);
+// //LagSpace = new ScalarXFEMToVectorFS(_tag,ScalarXFEMToVectorFS::VECTOR_X,ScalarXFEMToVectorFS::VECTOR_Y,_TagEnrichedVertex,_funcEnrichment);
+//
+// // we first do all fixations. the behavior of the dofManager is to
+// // give priority to fixations : when a dof is fixed, it cannot be
+// // numbered afterwards
+//
+//
+//
+//// for (int i = 0; i < elasticFields.size(); ++i){
+//// groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
+//// for ( ; it != elasticFields[i].g->end(); ++it)
+//// {
+//// MElement *e = *it;
+//// std::cout<<"elasticfield e num : "<<e->getNum()<<"\n";
+//// if (e->getParent()) std::cout<<"e->getParent : "<<e->getParent()->getNum()<<"\n";
+//// std::cout<<"vertex num 1: "<<e->getVertex(0)->getNum()<<"\n";
+//// std::cout<<"vertex num 2: "<<e->getVertex(1)->getNum()<<"\n";
+//// std::cout<<"vertex num 3: "<<e->getVertex(2)->getNum()<<"\n";
+//// }
+//// }
+//
+//
+//
+//
+// std::cout << "Dirichlet BC"<< std::endl;
+// for (unsigned int i = 0; i < allDirichlet.size(); i++)
+// {
+//// std::set<MElement*>::const_iterator it = allDirichlet[i].g->begin();
+//// for (;it!=allDirichlet[i].g->end();it++)
+//// {
+//// MElement *e = (*it);
+////// std::cout<<"dirch tag: "<<allDirichlet[i]._tag<<"\n";
+//////
+////// std::cout << "test" << std::endl;
+////
+//// std::cout<<"dirich e num : "<<e->getNum()<<"\n";
+//// std::cout<<"y pos : "<< e->getVertex(0)->y() <<"\n";
+//// e->getVertex(0)->y();
+////
+//// //std::cout<<"y pos : "<< e->getNum() <<"\n";
+//// //printf("y pos : %d",e->getVertex(1)->y());
+//// }
+//
+// FilterDofComponent filter(allDirichlet[i]._comp);
+// FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,allDirichlet[i]._f,filter);
+// }
+//
+// // we number the dofs : when a dof is numbered, it cannot be numbered
+// // again with another number.
+// for (unsigned int i = 0; i < elasticFields.size(); ++i)
+// {
+// NumberDofs(*LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),*pAssembler);
+// }
+//
+// // Now we start the assembly process
+// // First build the force vector
+//
+// GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
+// std::cout << "Neumann BC"<< std::endl;
+// for (unsigned int i = 0; i < allNeumann.size(); i++)
+// {
+// LoadTerm<SVector3> Lterm(*LagSpace,allNeumann[i]._f);
+// Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler);
+// }
+//
+//// bulk material law
+//
+// GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
+// for (unsigned int i = 0; i < elasticFields.size(); i++)
+// {
+// IsotropicElasticTerm Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu);
+//// LaplaceTerm<SVector3,SVector3> Eterm(*LagSpace);
+// Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
+// }
+//
+// printf("-- done assembling!\n");
+// std::cout<<"Dof Number : " << pAssembler->sizeOfR() <<"\n";
+// lsys->systemSolve();
+// printf("-- done solving!\n");
+//
+// double energ=0;
+// for (unsigned int i = 0; i < elasticFields.size(); i++)
+// {
+// SolverField<SVector3> Field(pAssembler, LagSpace);
+// IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
+// BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
+// Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
+// }
+// printf("elastic energy=%f\n",energ);
+// for (int i=0;i<pAssembler->sizeOfR();i++) std::cout<< lsys->getFromSolution(i) << "\n";
+
+}
+
+
+
+PView* OctreeSolver::buildDisplacementView (const std::string &postFileName)
+{
+
+ std::set<MVertex*> v;
+ for (unsigned int i = 0; i < elasticFields.size(); ++i)
+ {
+ for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
+ {
+ MElement *e=*it;
+ if (e->getParent()) e=e->getParent();
+ for (int j = 0; j < e->getNumVertices(); ++j) v.insert(e->getVertex(j));
+ }
+ }
+ std::map<int, std::vector<double> > data;
+ SolverField<SVector3> Field(pAssembler, LagSpace);
+ for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
+ {
+ SVector3 val;
+ MPoint p(*it);
+ Field.f(&p,0,0,0,val);
+ std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=val(2);
+// std::cout<<"ver num : " << (*it)->getNum() <<"\n" ;
+// std::cout<<"( "<<vec[0]<<" "<<vec[1]<<" "<<vec[2]<<" )\n";
+ data[(*it)->getNum()]=vec;
+ }
+ PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
+ return pv;
+
+}
+
diff --git a/contrib/arc/Classes/OctreeSolver.h b/contrib/arc/Classes/OctreeSolver.h
new file mode 100644
index 0000000000..4298f1e7c9
--- /dev/null
+++ b/contrib/arc/Classes/OctreeSolver.h
@@ -0,0 +1,53 @@
+//
+// Description :
+//
+//
+// Author: <Boris Sedji>, 01/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+
+#ifndef _OCTREE_SOLVER_H_
+#define _OCTREE_SOLVER_H_
+
+#include "elasticitySolver.h"
+#include "simpleFunction.h"
+//#include "QuadtreeLSImage.h"
+
+#include "gLevelSetMesh.cpp"
+#include "FuncGradDisc.h"
+
+class OctreeSolver : public elasticitySolver
+{
+ protected :
+
+ // map containing the tag of vertex and enriched status
+ std::set<int > _TagEnrichedVertex;
+ // enriched comp
+ std::set<int> _EnrichComp;
+ // simple multiplying function enrichment to enrich the space function
+ simpleFunction<double> *_funcEnrichment;
+ // level set value
+ std::map<int, double > _LevelSetValue;
+ // level set
+ gLevelSetMesh *_ls;
+
+ public :
+
+ OctreeSolver(int tag) : elasticitySolver(tag) {}
+ ~OctreeSolver();
+ //
+ virtual void readInputFile(const std::string &fn);
+ // create a GModel and determine de dimension of mesh in meshFileName
+ virtual void setMesh(const std::string &meshFileName);
+ // set the map _LevelSetValue
+ void setLevelSetValue(const std::string &LevelSetValueFileName);
+ // system solve, read the .dat file, fill tagEnrichedVertex, fill funcEnrichment, solve
+ virtual void solve();
+ virtual PView *buildDisplacementView(const std::string &postFileName);
+};
+
+#endif
+
diff --git a/contrib/arc/Classes/gLevelSetMesh.cpp b/contrib/arc/Classes/gLevelSetMesh.cpp
new file mode 100644
index 0000000000..c831283971
--- /dev/null
+++ b/contrib/arc/Classes/gLevelSetMesh.cpp
@@ -0,0 +1,21 @@
+//
+// Description :
+//
+//
+// Author: <Boris Sedji>, 01/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "gLevelSetMesh.h"
+
+#include "SPoint3.h"
+#include "MElement.h"
+#include "MVertex.h"
+#include "GModel.h"
+
+double gLevelSetMesh::getVertexValue(MVertex *v) const
+{
+ return (*_LevelSetValue)[v->getNum()];
+}
diff --git a/contrib/arc/Classes/gLevelSetMesh.h b/contrib/arc/Classes/gLevelSetMesh.h
new file mode 100644
index 0000000000..57ed49faf0
--- /dev/null
+++ b/contrib/arc/Classes/gLevelSetMesh.h
@@ -0,0 +1,61 @@
+//
+// Description :
+//
+//
+// Author: <Boris Sedji>, 01/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef _GLEVELSET_MESH_H_
+#define _GLEVELSET_MESH_H_
+
+#include "DILevelset.h"
+#include "MVertex.h"
+#include "GModel.h"
+
+class gLevelSetMesh : public gLevelsetPrimitive
+{
+ protected :
+
+ std::map<int, double > *_LevelSetValue;
+ GModel * _pModel;
+
+ public :
+ gLevelSetMesh(std::map<int, double > * LevelSetValue, GModel *p, int &tag) : gLevelsetPrimitive(tag)
+ {
+ _LevelSetValue = LevelSetValue;
+ _pModel = p;
+ };
+ // return negative value inward and positive value outward
+ virtual double operator() (const double &x, const double &y, const double &z) const
+ {
+ SPoint3 p(x,y,z);
+ MElement *e = _pModel->getMeshElementByCoord(p);
+ if (e->getParent()) e = e->getParent();
+ int numV = e->getNumVertices();
+ double xyz[3] = {x,y,z};
+ double uvw[3];
+ e->xyz2uvw(xyz,uvw);
+ double s[20];
+ e->getShapeFunctions(uvw[0],uvw[1],uvw[2],s);
+ double ls = 0;
+
+// std::cout<<"xyz : "<< x << " "<<y <<" "<<z<<"\n";
+ for (int i = 0;i<numV;i++)
+ {
+// std::cout<<"Vertex xyz :"<< e->getVertex(i)->x() << " "<<e->getVertex(i)->y() <<" "<<e->getVertex(i)->z()<<"\n";
+// std::cout<<"Vertex LS : "<< getVertexValue(e->getVertex(i)) <<"\n";
+ ls = ls + s[i]*getVertexValue(e->getVertex(i));
+ }
+// std::cout<<"=> LS : "<< ls <<"\n";
+ return ls;
+ }
+ virtual double getVertexValue(MVertex *v) const;
+ int type() const {return LSMESH;}
+};
+
+
+
+#endif
diff --git a/contrib/arc/Classes/xFemFunctionSpace.cpp b/contrib/arc/Classes/xFemFunctionSpace.cpp
new file mode 100644
index 0000000000..06fc20fadc
--- /dev/null
+++ b/contrib/arc/Classes/xFemFunctionSpace.cpp
@@ -0,0 +1,220 @@
+//
+// Description : Filtered and xFem function space definition
+//
+//
+// Author: <Eric Bechet>::<Boris Sedji>, 02/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "xFemFunctionSpace.h"
+
+
+template <class T> void xFemFS<T>::f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals)
+{
+ // We need parent parameters
+ MElement * elep;
+ if (ele->getParent()) elep = ele->getParent();
+ else elep = ele;
+
+ // Get the spacebase valsd
+ std::vector<ValType> valsd;
+ xFemFunctionSpace<T>::_spacebase->f(elep,u,v,w,valsd);
+
+ int nbdofs=valsd.size();
+ int curpos=vals.size(); // if in composite function space
+ vals.reserve(curpos+nbdofs);
+
+ // then enriched dofs so the order is ex:(a2x,a2y,a3x,a3y)
+ if (nbdofs>0) // if enriched
+ {
+ // Enrichment function calcul
+ SPoint3 p;
+ elep->pnt(u, v, w, p); // parametric to cartesian coordinates
+ double func;
+ func = (*_funcEnrichment)(p.x(), p.y(), p.z(),elep);
+ for (int i=0 ;i<nbdofs;i++)
+ {
+ vals.push_back(valsd[i]*func);
+ }
+ }
+}
+
+template <class T> void xFemFS<T>::gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+{
+
+ // We need parent parameters
+ MElement * elep;
+ if (ele->getParent()) elep = ele->getParent();
+ else elep = ele;
+
+ // Get the spacebase gradsd
+ std::vector<GradType> gradsd;
+ xFemFunctionSpace<T>::_spacebase->gradf(elep,u,v,w,gradsd);
+
+
+ int nbdofs=gradsd.size();
+
+ // We need spacebase valsd to compute total gradient
+ std::vector<ValType> valsd;
+ xFemFunctionSpace<T>::_spacebase->f(elep,u,v,w,valsd);
+
+ int curpos=grads.size(); // if in composite function space
+ grads.reserve(curpos+nbdofs);
+
+ // then enriched dofs so the order is ex:(a2x,a2y,a3x,a3y)
+ if (nbdofs>0) // if enriched
+ {
+ double df[3];
+ SPoint3 p;
+ elep->pnt(u, v, w, p);
+ _funcEnrichment->gradient (p.x(), p.y(),p.z(),df[0],df[1],df[2],elep);
+ ValType gradfuncenrich(df[0],df[1],df[2]);
+
+ // Enrichment function calcul
+
+ double func;
+ func = (*_funcEnrichment)(p.x(), p.y(), p.z(),elep);
+
+ for (int i=0 ;i<nbdofs;i++)
+ {
+ GradType GradFunc;
+ tensprod(valsd[i],gradfuncenrich,GradFunc);
+ grads.push_back(gradsd[i]*func + GradFunc);
+ }
+ }
+}
+
+template <class T> int xFemFS<T>::getNumKeys(MElement *ele)
+{
+ MElement * elep;
+ if (ele->getParent()) elep = ele->getParent();
+ else elep = ele;
+ int nbdofs = xFemFunctionSpace<T>::_spacebase->getNumKeys(ele);
+ return nbdofs;
+}
+
+template <class T> void xFemFS<T>::getKeys(MElement *ele, std::vector<Dof> &keys)
+{
+ MElement * elep;
+ if (ele->getParent()) elep = ele->getParent();
+ else elep = ele;
+
+ int normalk=xFemFunctionSpace<T>::_spacebase->getNumKeys(elep);
+
+ std::vector<Dof> bufk;
+ bufk.reserve(normalk);
+ xFemFunctionSpace<T>::_spacebase->getKeys(elep,bufk);
+ int normaldofs=bufk.size();
+ int nbdofs = normaldofs;
+
+ int curpos=keys.size();
+ keys.reserve(curpos+nbdofs);
+
+ // get keys so the order is ex:(a2x,a2y,a3x,a3y)
+ // enriched dof tagged with ( i2 -> i2 + 1 )
+ for (int i=0 ;i<nbdofs;i++)
+ {
+ int i1,i2;
+ Dof::getTwoIntsFromType(bufk[i].getType(), i1,i2);
+ keys.push_back(Dof(bufk[i].getEntity(),Dof::createTypeWithTwoInts(i1,i2+1)));
+ }
+}
+
+
+// Filtered function space
+//
+
+template <class T> void FilteredFS<T>::f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals)
+{
+ // We need parent parameters
+ MElement * elep;
+ if (ele->getParent()) elep = ele->getParent();
+ else elep = ele;
+
+ std::vector<ValType> valsd;
+
+ FilteredFunctionSpace<T,FilterNodeEnriched>::_spacebase->f(elep,u,v,w,valsd);
+
+ int normalk=FilteredFunctionSpace<T,FilterNodeEnriched>::_spacebase->getNumKeys(elep);
+ std::vector<Dof> bufk;
+ bufk.reserve(normalk);
+ FilteredFunctionSpace<T,FilterNodeEnriched>::_spacebase->getKeys(elep,bufk);
+
+ for (int i=0;i<bufk.size();i++)
+ {
+ if ((*FilteredFunctionSpace<T,FilterNodeEnriched>::_filter)(bufk[i]))
+ vals.push_back(valsd[i]);
+ }
+
+}
+
+
+template <class T> void FilteredFS<T>::gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+{
+ // We need parent parameters
+ MElement * elep;
+ if (ele->getParent()) elep = ele->getParent();
+ else elep = ele;
+
+ // Get space base gradsd
+ std::vector<GradType> gradsd;
+ FilteredFunctionSpace<T,FilterNodeEnriched>::_spacebase->gradf(elep,u,v,w,gradsd);
+
+ // Get numkeys
+ int normalk=FilteredFunctionSpace<T,FilterNodeEnriched>::_spacebase->getNumKeys(elep);
+ std::vector<Dof> bufk;
+ bufk.reserve(normalk);
+ FilteredFunctionSpace<T,FilterNodeEnriched>::_spacebase->getKeys(elep,bufk);
+
+ for (int i=0;i<bufk.size();i++)
+ {
+ if ( (*FilteredFunctionSpace<T,FilterNodeEnriched>::_filter)(bufk[i]))
+ grads.push_back(gradsd[i]);
+ }
+
+}
+
+template <class T> int FilteredFS<T>::getNumKeys(MElement *ele)
+{
+ MElement *elep;
+ if (ele->getParent()) elep = ele->getParent();
+ else elep = ele;
+
+ int nbdofs = 0;
+
+ int normalk=FilteredFunctionSpace<T,FilterNodeEnriched>::_spacebase->getNumKeys(elep);
+ std::vector<Dof> bufk;
+ bufk.reserve(normalk);
+ FilteredFunctionSpace<T,FilterNodeEnriched>::_spacebase->getKeys(elep,bufk);
+
+ for (int i=0;i<bufk.size();i++)
+ {
+ if ( (*FilteredFunctionSpace<T,FilterNodeEnriched>::_filter)(bufk[i]))
+ nbdofs = nbdofs + 1;
+ }
+ return nbdofs;
+}
+
+template <class T> void FilteredFS<T>::getKeys(MElement *ele, std::vector<Dof> &keys)
+{
+ MElement * elep;
+ if (ele->getParent()) elep = ele->getParent();
+ else elep = ele;
+
+ int normalk=FilteredFunctionSpace<T,FilterNodeEnriched>::_spacebase->getNumKeys(elep);
+
+ std::vector<Dof> bufk;
+ bufk.reserve(normalk);
+ FilteredFunctionSpace<T,FilterNodeEnriched>::_spacebase->getKeys(elep,bufk);
+ int normaldofs=bufk.size();
+ int nbdofs = 0;
+
+ for (int i=0;i<bufk.size();i++)
+ {
+ if ( (*FilteredFunctionSpace<T,FilterNodeEnriched>::_filter)(bufk[i]))
+ keys.push_back(bufk[i]);
+ }
+}
+
diff --git a/contrib/arc/Classes/xFemFunctionSpace.h b/contrib/arc/Classes/xFemFunctionSpace.h
new file mode 100644
index 0000000000..b081870d5a
--- /dev/null
+++ b/contrib/arc/Classes/xFemFunctionSpace.h
@@ -0,0 +1,111 @@
+//
+// Description : Filtered and xFem function space definition
+//
+//
+// Author: <Eric Bechet>::<Boris Sedji>, 02/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef _XFEMFS_H_
+#define _XFEMFS_H_
+
+#include "simpleFunction.h"
+
+class FilterNodeEnriched
+{
+
+ private :
+
+ std::set<int> *_TagEnrichedVertex;
+ std::set<int> * _EnrichComp;
+
+ public :
+
+ FilterNodeEnriched(std::set<int > * TagEnrichedVertex , std::set<int> * EnrichComp)
+ {
+ _TagEnrichedVertex = TagEnrichedVertex;
+ _EnrichComp = EnrichComp;
+ }
+
+ virtual bool operator () (Dof & key) const
+ {
+ std::set<int>::iterator it1;
+ std::set<int>::iterator it2;
+ int i1,i2;
+ Dof::getTwoIntsFromType(key.getType(), i1,i2);
+ it2 = _EnrichComp->find(i1);
+ it1 = _TagEnrichedVertex->find(key.getEntity());
+ if (it1!=_TagEnrichedVertex->end() & it2 != _EnrichComp->end())
+ {
+ return true;
+ }
+ else return false;
+ }
+
+ //std::vector<int> * getEnrichComp(){return _EnrichComp;}
+
+// void SetEnrichedVertex(MElement *elep, std::vector<int> & EnrichedVertex,int &nbdofs)
+// {
+// EnrichedVertex.clear();
+// nbdofs = 0;
+// for (int i=0 ;i<elep->getNumVertices();i++)
+// {
+// std::set<int>::iterator it;
+// it = _TagEnrichedVertex->find(elep->getVertex(i)->getNum());
+// if (it!=_TagEnrichedVertex->end())
+// {
+// EnrichedVertex.push_back(i);
+// nbdofs = nbdofs + 1*_EnrichComp->size(); // enriched dof
+// }
+// }
+// }
+};
+
+
+template<class T>
+class xFemFS : public xFemFunctionSpace<T>
+{
+
+ public :
+
+ typedef typename TensorialTraits<T>::ValType ValType;
+ typedef typename TensorialTraits<T>::GradType GradType;
+
+ protected:
+
+ simpleFunction<double> *_funcEnrichment;
+
+ public:
+
+ xFemFS(FunctionSpace<T>* spacebase, simpleFunction<double> *funcEnrichment) : xFemFunctionSpace<T>(spacebase),_funcEnrichment(funcEnrichment)
+ {};
+
+ virtual void f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
+ virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
+ virtual int getNumKeys(MElement *ele);
+ virtual void getKeys(MElement *ele, std::vector<Dof> &keys);
+};
+
+
+
+template<class T>
+class FilteredFS : public FilteredFunctionSpace<T,FilterNodeEnriched>
+{
+
+ public :
+
+ typedef typename TensorialTraits<T>::ValType ValType;
+ typedef typename TensorialTraits<T>::GradType GradType;
+
+ public:
+ FilteredFS(FunctionSpace<T>* spacebase, FilterNodeEnriched *filter) : FilteredFunctionSpace<T,FilterNodeEnriched> (spacebase,filter)
+ {}
+ virtual void f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
+ virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
+ virtual int getNumKeys(MElement *ele);
+ virtual void getKeys(MElement *ele, std::vector<Dof> &keys);
+};
+
+#endif
--
GitLab