diff --git a/contrib/arc/Classes/FuncHeaviside.h b/contrib/arc/Classes/FuncHeaviside.h
new file mode 100644
index 0000000000000000000000000000000000000000..4cd9895d362bd42f074482ef61660a52de2dfdfc
--- /dev/null
+++ b/contrib/arc/Classes/FuncHeaviside.h
@@ -0,0 +1,39 @@
+//
+// Description : Heaviside function based on level set discontinuity description
+//
+//
+// Author: <Boris Sedji>, 01/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef _FUNCHEAVISIDE_H_
+#define _FUNCHEAVISIDE_H_
+
+#include "simpleFunction.h"
+#include "DILevelset.h"
+
+
+class FuncHeaviside : public simpleFunction<double> {
+ private :
+
+ gLevelset *_ls;
+
+ public :
+ FuncHeaviside(gLevelset *ls){
+ _ls = ls;
+ }
+ virtual double operator () (double x,double y,double z) const {
+ if (_ls->isInsideDomain (x,y,z)){
+ return 1 ;
+ }else{
+ return -1;
+ }
+ }
+ virtual void gradient (double x, double y, double z,
+ double & dfdx, double & dfdy, double & dfdz) const
+ { dfdx = dfdy = dfdz = 0.0; }
+};
+
+#endif
diff --git a/contrib/arc/Classes/VectorXFEMFS.cpp b/contrib/arc/Classes/VectorXFEMFS.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0257484186ad2be7567a123cc0b61286019aa56e
--- /dev/null
+++ b/contrib/arc/Classes/VectorXFEMFS.cpp
@@ -0,0 +1,177 @@
+//
+// Description : XFEM elasticity solver implementation
+//
+//
+// Author: <Eric Bechet>::<Boris Sedji>, 01/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "VectorXFEMFS.h"
+#include "SPoint3.h"
+
+void ScalarLagrangeToXfemFS::f(MVertex *ver, std::vector<ValType> &vals)
+{
+ if (!(*_TagEnrichedVertex)[ver->getNum()]) vals.push_back(1.0);
+ else
+ {
+ double func = (*_funcEnrichment)(ver->x(),ver->y(),ver->z());
+ vals.push_back(1.0);
+ vals.push_back(func);
+ }
+}
+
+void ScalarLagrangeToXfemFS::f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
+{
+ // Enrichment function calcul
+ SPoint3 p;
+ ele->pnt(u, v, w, p);
+ double func;
+ func = (*_funcEnrichment)(p.x(), p.y(), p.z());
+
+ // We need parent parameters
+ if (ele->getParent()) ele = ele->getParent();
+
+ int ndofs,normaldofs;
+ ndofs = ele->getNumVertices();
+ normaldofs = ndofs;
+
+ for (int i=0 ;i<ele->getNumVertices();i++)
+ {
+ if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()]) ndofs = ndofs + 1; // enriched dof
+ }
+
+ int curpos=vals.size();
+ vals.resize(curpos+ndofs);
+
+ // normal shape functions
+ ele->getShapeFunctions(u, v, w, &(vals[curpos]));
+
+ int k=0;
+ for (int i=0 ;i<ele->getNumVertices();i++)
+ {
+ if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()])
+ {
+ vals[curpos+normaldofs+k] = vals[curpos+i]*func; // enriched dof
+ k++;
+ }
+ }
+
+};
+
+void ScalarLagrangeToXfemFS::gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
+{
+ // Enrichment function calcul
+ SPoint3 p;
+ ele->pnt(u, v, w, p);
+ double func;
+ func = (*_funcEnrichment)(p.x(), p.y(), p.z());
+
+ // We need parent parameters
+ if (ele->getParent()) ele = ele->getParent();
+
+ int ndofs,normaldofs;
+ ndofs= ele->getNumVertices();
+ normaldofs = ndofs; // normal dofs
+
+ for (int i=0 ;i< (ele->getNumVertices());i++)
+ {
+ if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()]) ndofs = ndofs + 1; // enriched dof
+ }
+
+ int curpos = grads.size();
+ grads.reserve(curpos+ndofs);
+ double gradsuvw[256][3];
+ ele->getGradShapeFunctions(u, v, w, gradsuvw);
+
+ int k = 0;
+ for (int i=0 ;i<(ele->getNumVertices());i++)
+ {
+ if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()]) // enriched dof
+ {
+ gradsuvw[curpos+normaldofs+k][0] = gradsuvw[curpos+i][0]*func;
+ gradsuvw[curpos+normaldofs+k][1] = gradsuvw[curpos+i][1]*func;
+ gradsuvw[curpos+normaldofs+k][2] = gradsuvw[curpos+i][2]*func;
+ k++;
+ }
+ }
+
+ double jac[3][3];
+ double invjac[3][3];
+ const double detJ = ele->getJacobian(u, v, w, jac);
+ // redondant : on fait cet appel a l'exterieur
+ inv3x3(jac, invjac);
+ for (int i=0;i<ndofs;++i)
+ grads.push_back(GradType(
+ invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2],
+ invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2],
+ invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2]
+ ));
+
+
+};
+
+int ScalarLagrangeToXfemFS::getNumKeys(MElement *ele)
+{
+ // We need parent parameters
+ if (ele->getParent()) ele = ele->getParent();
+ int ndofs;
+ ndofs = ele->getNumVertices();
+ {
+ for (int i=0 ;i<(ele->getNumVertices());i++)
+ {
+ if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()]) ndofs = ndofs + 1;
+ }
+ }
+ return ndofs;
+}
+
+int ScalarLagrangeToXfemFS::getNumKeys(MVertex *ver)
+{
+ if (!(*_TagEnrichedVertex)[ver->getNum()]) return 1 ;
+ else return 2; // if enriched vertex, there is two dof at this vertex (one dimension)
+}
+
+void ScalarLagrangeToXfemFS::getKeys(MElement *ele, std::vector<Dof> &keys)
+{
+ if (ele->getParent()) ele = ele->getParent();
+ // vector of additionals keys
+ std::vector<Dof> SupKeys;
+ // vector of all vertex dof keys
+ std::vector<Dof> VertexKeys;
+
+ int ndofs,normaldofs;
+ ndofs = ele->getNumVertices();
+
+ for (int i=0 ;i<(ele->getNumVertices());i++)
+ {
+ if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()]) ndofs = ndofs + 1;
+ }
+
+ keys.reserve(keys.size()+ndofs);
+
+ // the order of dof in one element of three vertex is u1x u2x u3x a1x a2x a3x u1y, etc...
+ for (int i=0;i< (ele->getNumVertices()) ;++i)
+ {
+ getKeys(ele->getVertex(i), VertexKeys);
+ if (VertexKeys.size()>1) SupKeys.push_back(VertexKeys[1]);
+ keys.push_back(VertexKeys[0]);
+ VertexKeys.clear();
+ }
+
+ for (int i=0;i< SupKeys.size() ;++i) keys.push_back(SupKeys[i]);
+
+}
+
+
+void ScalarLagrangeToXfemFS::getKeys(MVertex *ver, std::vector<Dof> &keys)
+{
+ if ((*_TagEnrichedVertex)[ver->getNum()])
+ {
+ keys.push_back(Dof(ver->getNum(), _iField));
+ keys.push_back(Dof(ver->getNum(), _iField+1)); // we tag the additional dof with number field+1
+ }
+ else keys.push_back(Dof(ver->getNum(), _iField));
+};
+
diff --git a/contrib/arc/Classes/VectorXFEMFS.h b/contrib/arc/Classes/VectorXFEMFS.h
new file mode 100644
index 0000000000000000000000000000000000000000..92132fef335dcef94572e90484aa7284147869aa
--- /dev/null
+++ b/contrib/arc/Classes/VectorXFEMFS.h
@@ -0,0 +1,97 @@
+//
+// Description : From Scalar Lagrange Function Space To XFEM on Tagged Vertex
+//
+//
+// Author: <Eric Bechet>::<Boris Sedji>, 01/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "SVector3.h"
+#include "STensor3.h"
+#include <vector>
+#include <iterator>
+#include <iostream>
+#include "Numeric.h"
+#include "MElement.h"
+#include "dofManager.h"
+
+
+#include "functionSpace.h"
+#include "simpleFunction.h"
+
+
+class ScalarLagrangeToXfemFS : public ScalarLagrangeFunctionSpace{
+ public:
+
+ typedef TensorialTraits<double>::ValType ValType;
+ typedef TensorialTraits<double>::GradType GradType;
+ typedef TensorialTraits<double>::HessType HessType;
+ typedef TensorialTraits<double>::DivType DivType;
+ typedef TensorialTraits<double>::CurlType CurlType;
+
+ protected:
+
+ std::map<int , bool> *_TagEnrichedVertex;
+ simpleFunction<double> *_funcEnrichment;
+
+ public:
+//
+ ScalarLagrangeToXfemFS(int i, std::map<int , bool> &TagEnrichedVertex, simpleFunction<double> *funcEnrichment) : ScalarLagrangeFunctionSpace(i)
+ {
+ _TagEnrichedVertex = & TagEnrichedVertex;
+ _funcEnrichment = funcEnrichment;
+ }
+
+ virtual int getId(void) const {return ScalarLagrangeFunctionSpace::_iField;};
+ // Shape function value in element at uvw
+ virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals);
+ // Shape function value at vertex
+ virtual void f(MVertex *ver, std::vector<ValType> &vals) ;
+ // Grad Shape function value in element at uvw
+ virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) ;
+ // Shape function Number for element (in one dimension (scalar))
+ virtual int getNumKeys(MElement *ele);
+ // Dof keys for the element
+ virtual void getKeys(MElement *ele, std::vector<Dof> &keys) ;
+ // Shape function Number associate with the vertex
+ virtual int getNumKeys(MVertex *ver);
+ // Get Dof keys for the vertex (in one dimension (scalar))
+ virtual void getKeys(MVertex *ver, std::vector<Dof> &keys);
+};
+
+
+class ScalarXFEMToVectorFS : public ScalarToAnyFunctionSpace<SVector3>
+{
+ protected:
+
+ static const SVector3 BasisVectors[3];
+
+ public:
+ enum Along { VECTOR_X=0, VECTOR_Y=1, VECTOR_Z=2 };
+
+ typedef TensorialTraits<SVector3>::ValType ValType;
+ typedef TensorialTraits<SVector3>::GradType GradType;
+ typedef TensorialTraits<SVector3>::HessType HessType;
+ typedef TensorialTraits<SVector3>::DivType DivType;
+ typedef TensorialTraits<SVector3>::CurlType CurlType;
+
+ ScalarXFEMToVectorFS(int id , std::map<int , bool> & TagEnrichedVertex , simpleFunction<double> * funcEnrichment) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment),
+ SVector3(1.,0.,0.),VECTOR_X, SVector3(0.,1.,0.),VECTOR_Y,SVector3(0.,0.,1.),VECTOR_Z) {}
+
+ ScalarXFEMToVectorFS(int id,Along comp1, std::map<int , bool> &TagEnrichedVertex , simpleFunction<double> * funcEnrichment) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment),
+ BasisVectors[comp1],comp1) {}
+
+ ScalarXFEMToVectorFS(int id,Along comp1,Along comp2, std::map<int , bool> &TagEnrichedVertex,simpleFunction<double> *funcEnrichment) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment),
+ BasisVectors[comp1],comp1, BasisVectors[comp2],comp2) {}
+
+ ScalarXFEMToVectorFS(int id,Along comp1,Along comp2, Along comp3, std::map<int , bool> &TagEnrichedVertex,simpleFunction<double> *funcEnrichment) :
+ ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment),
+ BasisVectors[comp1],comp1, BasisVectors[comp2],comp2, BasisVectors[comp3],comp3) {}
+};
+
+const SVector3 ScalarXFEMToVectorFS::BasisVectors[3]={SVector3(1,0,0),SVector3(0,1,0),SVector3(0,0,1)};
diff --git a/contrib/arc/Classes/XFEMelasticitySolver.cpp b/contrib/arc/Classes/XFEMelasticitySolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e2624c411080387438bfddcd091a0a5252633060
--- /dev/null
+++ b/contrib/arc/Classes/XFEMelasticitySolver.cpp
@@ -0,0 +1,183 @@
+//
+// Description : XFEM elasticity solver, element space function enriched on tagged vertex
+//
+//
+// Author: <Eric Bechet>::<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"
+#endif
+
+
+#include "DILevelset.h"
+#include "VectorXFEMFS.cpp"
+#include "XFEMelasticitySolver.h"
+#include "FuncHeaviside.h"
+
+
+XFEMelasticitySolver::~XFEMelasticitySolver()
+{
+ delete _funcEnrichment;
+}
+
+
+void XFEMelasticitySolver::setMesh(const std::string &meshFileName)
+{
+ pModel = new GModel();
+ pModel->readMSH(meshFileName.c_str());
+ _dim = pModel->getNumRegions() ? 3 : 2;
+}
+
+void XFEMelasticitySolver::solve(){
+
+#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){
+ groupOfElements::elementContainer::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[e->getParent()->getVertex(k)->getNum()] = true;
+ }
+ }
+ else
+ {
+ for (int k = 0; k < e->getNumVertices(); ++k)
+ {
+ _TagEnrichedVertex[e->getVertex(k)->getNum()] = false;
+ }
+ }
+ }
+ }
+
+ // 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;
+ if (_dim==3) LagSpace=new ScalarXFEMToVectorFS(_tag,_TagEnrichedVertex,_funcEnrichment);
+ if (_dim==2) 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 (unsigned int i = 0; i < allDirichlet.size(); i++)
+ {
+ FilterDofComponent filter(allDirichlet[i]._comp);
+ if (allDirichlet[i].onWhat==BoundaryCondition::ON_VERTEX)
+ FixNodalDofs(*LagSpace,allDirichlet[i].g->vbegin(),allDirichlet[i].g->vend(),*pAssembler,allDirichlet[i]._f,filter);
+ else
+ 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);
+
+ for (unsigned int i = 0; i < allNeumann.size(); i++)
+ {
+ LoadTerm<SVector3> Lterm(*LagSpace,allNeumann[i]._f);
+ if (allNeumann[i].onWhat==BoundaryCondition::ON_VERTEX)
+ Assemble(Lterm,*LagSpace,allNeumann[i].g->vbegin(),allNeumann[i].g->vend(),*pAssembler);
+ else
+ 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);
+ 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);
+}
+
+
+
+PView* XFEMelasticitySolver::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;
+ Field.f(*it,val);
+ std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=0;
+ data[(*it)->getNum()]=vec;
+ }
+ PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
+ return pv;
+}
diff --git a/contrib/arc/Classes/XFEMelasticitySolver.h b/contrib/arc/Classes/XFEMelasticitySolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..2c587069c60f619cb09f3369c55e788f18afb6b9
--- /dev/null
+++ b/contrib/arc/Classes/XFEMelasticitySolver.h
@@ -0,0 +1,36 @@
+//
+// Description : XFEM elasticity solver, element space function enriched on tagged vertex
+//
+//
+// Author: <Eric Bechet>::<Boris Sedji>, 01/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#ifndef _XFEMELASTICITY_SOLVER_H_
+#define _XFEMELASTICITY_SOLVER_H_
+
+#include "elasticitySolver.h"
+#include "simpleFunction.h"
+
+class XFEMelasticitySolver : public elasticitySolver
+{
+ protected :
+ // map containing the tag of vertex and enriched status
+ std::map<int , bool> _TagEnrichedVertex;
+ // simple multiplying function enrichment to enrich the space function
+ simpleFunction <double> *_funcEnrichment;
+
+ public :
+
+ XFEMelasticitySolver(int tag) : elasticitySolver(tag) {}
+ ~XFEMelasticitySolver();
+ // create a GModel and determine de dimension of mesh in meshFileName
+ virtual void setMesh(const std::string &meshFileName);
+ // 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/ImportLS2dImage.cpp b/contrib/arc/ImportLS2dImage.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..34cc9e82f6daaf25c108e7af5687672471ed545d
--- /dev/null
+++ b/contrib/arc/ImportLS2dImage.cpp
@@ -0,0 +1,84 @@
+//
+// Description : Import Float Level Set Image in Gmsh QuadTree Mesh
+//
+//
+// Author: <Boris Sedji>, 12/2009
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+
+
+#include "itkImage.h"
+#include "itkImageFileReader.h"
+#include "itkImageFileWriter.h"
+#include <stdio.h>
+#include "GModel.h"
+#include "Classes/QuadtreeLSImage.cpp"
+
+
+
+int main( int argc, char *argv[] )
+{
+ if( argc < 4 )
+ {
+ std::cerr << "Missing Parameters " << std::endl;
+ std::cerr << "Usage: " << argv[0];
+ std::cerr << " inputLevelSet2dImage MeshSizeMax MeshSizeMin"<< std::endl;
+ return 1;
+ }
+
+ typedef float PixelTypeFloat;
+ const unsigned int Dimension = 2;
+ typedef itk::Image< PixelTypeFloat, Dimension > ImageTypeFloat;
+ typedef itk::ImageFileReader< ImageTypeFloat > ImageReaderTypeFloat;
+ ImageReaderTypeFloat::Pointer reader = ImageReaderTypeFloat::New();
+
+
+ reader->SetFileName(argv[1]);
+
+ ImageTypeFloat::Pointer image = reader->GetOutput();
+ image->Update();
+
+ ImageTypeFloat::RegionType region;
+ region = image->GetLargestPossibleRegion ();
+
+ std::cout<<"\nImage dimensions : " << region.GetSize(0) << " x " << region.GetSize(1);
+
+ QuadtreeLSImage quadtree(image);
+
+ int sizemax = atoi(argv[2]);
+ int sizemin = atoi(argv[3]);
+
+ // Mesh with conditions on mesh size
+ quadtree.Mesh(sizemax,sizemin);
+
+ // Smoothing mesh to avoid too much generation between adjacent elements
+ bool statut=false;
+ int k = 0;
+ std::cout<<"\nLeaf Number : "<< (quadtree.GetLeafNumber())<<"\n";
+ while((!statut) & (k < 20)){
+ statut = quadtree.Smooth();
+ std::cout<<"\nk : "<< k;
+ std::cout<<"\nsmoothed : " << (int)statut;
+ std::cout<<"\nLeaf Number : "<< (quadtree.GetLeafNumber())<<"\n";
+ k++;
+ }
+
+ // Create GModel with the octree mesh
+ GModel *m = quadtree.CreateGModel();
+
+ // Write a .msh file with the mesh
+ std::string ModelName = "QuadtreeMesh.msh" ;
+ m->writeMSH(ModelName,2.1,false,false);
+
+ // Write a .msh file with the level set values as postview data
+ PView *pv = quadtree.CreateLSPView(m);
+ bool useadapt = true;
+ pv->getData(useadapt)->writeMSH("LSPView.msh", false);
+
+ std::cout<<"\n";
+
+ return 0;
+}
diff --git a/contrib/arc/XFEMElasticity.cpp b/contrib/arc/XFEMElasticity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..17b5246974172e11225c7f3adf8f16ac56f2cf88
--- /dev/null
+++ b/contrib/arc/XFEMElasticity.cpp
@@ -0,0 +1,55 @@
+//
+// Description : XFEM elasticity solver main
+//
+//
+// Author: <Eric Bechet>::<Boris Sedji>, 01/2010
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+#include "Gmsh.h"
+#include "elasticitySolver.h"
+#include "PView.h"
+#include "PViewData.h"
+#include "highlevel.h"
+#include "groupOfElements.h"
+#include <iterator>
+
+#include "Classes/XFEMelasticitySolver.cpp"
+
+
+
+int main (int argc, char* argv[])
+{
+
+ if (argc < 2)
+ {
+ printf("usage : elasticity input_file_name yeahhh \n");
+ return -1;
+ }
+
+ GmshInitialize(argc, argv);
+ // globals are still present in Gmsh
+
+ // instanciate a solver
+ XFEMelasticitySolver mySolver (1000);
+
+ // read some input file
+ mySolver.readInputFile(argv[1]);
+
+ // solve the problem
+ mySolver.solve();
+
+ PView *pv = mySolver.buildDisplacementView("displacement");
+ pv->getData()->writeMSH("disp.msh", false);
+ delete pv;
+ pv = mySolver.buildElasticEnergyView("elastic energy");
+ pv->getData()->writeMSH("energ.msh", false);
+ delete pv;
+
+ // stop gmsh
+ GmshFinalize();
+
+}
+