diff --git a/Numeric/simpleFunction.h b/Numeric/simpleFunction.h
index 6cadc0eb27b007a9b74a9bc14b73abc57cbef424..4960b290d05a7e7a0d40e1a0cf28a3bed88572d1 100644
--- a/Numeric/simpleFunction.h
+++ b/Numeric/simpleFunction.h
@@ -6,6 +6,8 @@
#ifndef _SIMPLE_FUNCTION_H_
#define _SIMPLE_FUNCTION_H_
+#include "MElement.h"
+
template <class scalar>
class simpleFunction {
scalar _val;
@@ -14,8 +16,32 @@ class simpleFunction {
virtual ~simpleFunction(){}
virtual scalar operator () (double x, double y, double z) const { return _val; }
virtual void gradient (double x, double y, double z,
- scalar & dfdx, scalar & dfdy, scalar & dfdz) const
+ scalar & dfdx, scalar & dfdy, scalar & dfdz) const
{ dfdx = dfdy = dfdz = 0.0; }
+/* virtual scalar operator () (double x, double y, double z,MElement *e) const { return _val; }
+ virtual void gradient (double x, double y, double z,
+ scalar & dfdx, scalar & dfdy, scalar & dfdz, MElement *e) const
+ { dfdx = dfdy = dfdz = 0.0; }*/
+};
+
+template <class scalar>
+class simpleFunctionOnElement : public simpleFunction<scalar>
+{
+ MElement *_e;
+ public :
+ simpleFunctionOnElement(scalar val=0) : simpleFunction<scalar>(val),_e(0) {}
+ virtual ~simpleFunctionOnElement(){}
+ void setElement(MElement *e) {_e=e;}
+ MElement * getElement(void) const {return _e;}
+ MElement * getElement(double x, double y, double z) const
+ {
+ if (_e) return _e;
+ else
+ {// search
+ }
+ }
};
+
+
#endif
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 31a1fb782bc3b59d5b667a13840a382271b73a63..942a375092395e1364b4ce628c27d497ad06661f 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -14,6 +14,8 @@
#include "linearSystem.h"
#include "fullMatrix.h"
+#include <iostream>
+
class Dof{
private:
// v(x) = \sum_f \sum_i v_{fi} s^f_i(x)
@@ -49,8 +51,8 @@ template<class T> struct dofTraits
typedef T VecType;
typedef T MatType;
inline static void mult (VecType &r, const MatType &m, const VecType &v) { r = m*v;}
- inline static void neg (VecType &r) { r = -r;}
- inline static void setToZero (VecType &r) { r = 0;}
+ inline static void neg (VecType &r) { r = -r;}
+ inline static void setToZero (VecType &r) { r = 0;}
};
template<class T> struct dofTraits<fullMatrix<T> >
@@ -58,8 +60,8 @@ template<class T> struct dofTraits<fullMatrix<T> >
typedef fullMatrix<T> VecType;
typedef fullMatrix<T> MatType;
inline static void mult (VecType &r, const MatType &m, const VecType &v) { m.mult(v, r);}
- inline static void neg (VecType &r) { r.scale(-1.);}
- inline static void setToZero (VecType &r) { r.scale(0);}
+ inline static void neg (VecType &r) { r.scale(-1.);}
+ inline static void setToZero (VecType &r) { r.scale(0);}
};
/*
template<> struct dofTraits<fullVector<std::complex<double> > >
@@ -86,12 +88,13 @@ class dofManager{
private:
// general affine constraint on sub-blocks, treated by adding
// equations:
- // dataMat * DofVec = \sum_i dataMat_i * DofVec_i + dataVec
- std::map<std::pair<Dof, dataMat>, DofAffineConstraint<dataVec> > constraints;
+
+ // Dof = \sum_i dataMat_i x Dof_i + dataVec
+ std::map<Dof, DofAffineConstraint< dataVec > > constraints;
// fixations on full blocks, treated by eliminating equations:
// DofVec = dataVec
- std::map<Dof, dataVec> fixed ;
+ std::map<Dof, dataVec> fixed;
// initial conditions
std::map<Dof, std::vector<dataVec> > initial;
@@ -107,14 +110,15 @@ class dofManager{
linearSystem<dataMat> *_current;
public:
- dofManager() : _current(0) { }
dofManager(linearSystem<dataMat> *l) : _current(l) { _linearSystems["A"] = l; }
- dofManager(linearSystem<dataMat> *l1, linearSystem<dataMat> *l2) : _current(l1) {
+ dofManager(linearSystem<dataMat> *l1, linearSystem<dataMat> *l2) : _current(l1) {
_linearSystems.insert(std::make_pair("A", l1));
_linearSystems.insert(std::make_pair("B", l2));
+ //_linearSystems.["A"] = l1;
+ //_linearSystems["B"] = l2;
}
inline void fixDof(Dof key, const dataVec &value)
- {
+ {
fixed[key] = value;
}
inline void fixDof(long int ent, int type, const dataVec &value)
@@ -124,11 +128,18 @@ class dofManager{
inline void fixVertex(MVertex*v, int iComp, int iField, const dataVec &value)
{
fixDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField), value);
- }
+ }
inline void numberDof(Dof key)
{
- if(fixed.find(key) != fixed.end()) return;
- // if (constraints.find(key) != constraints.end()) return;
+ if(fixed.find(key) != fixed.end())
+ {
+ return;
+ }
+ // constraints
+ if (constraints.find(key) != constraints.end()){
+ return;
+ }
+ //
std::map<Dof, int> :: iterator it = unknown.find(key);
if (it == unknown.end()){
unsigned int size = unknown.size();
@@ -162,13 +173,27 @@ class dofManager{
if (it != unknown.end())
return _current->getFromSolution(it->second);
}
+ {
+ typename std::map<Dof, DofAffineConstraint< dataVec > >::const_iterator it = constraints.find(key);
+ if (it != constraints.end())
+ {
+ dataVec somme(0.0);
+ for (int i=0;i<(it->second).linear.size();i++)
+ {
+ std::map<Dof, int>::const_iterator itu = unknown.find(((it->second).linear[i]).first);
+ somme = somme + getDofValue(((it->second).linear[i]).first)*((it->second).linear[i]).second;
+ }
+ somme = somme + (it->second).shift;
+ return somme;
+ }
+ }
return dataVec(0.0);
}
inline void getDofValue(dataVec &v, int ent, int type) const
{
Dof key(ent, type);
- {
+ {
typename std::map<Dof, dataVec>::const_iterator it = fixed.find(key);
if (it != fixed.end()){
v = it->second;
@@ -182,8 +207,24 @@ class dofManager{
return;
}
}
+ {
+ typename std::map<Dof, DofAffineConstraint< dataVec > >::const_iterator it = constraints.find(key);
+ if (it != constraints.end())
+ {
+ dofTraits<T>::setToZero(v);
+ for (int i=0;i<(it->second).linear.size();i++)
+ {
+ std::map<Dof, int>::const_iterator itu = unknown.find(((it->second).linear[i]).first);
+ v = v + getDofValue(((it->second).linear[i]).first)*((it->second).linear[i]).second;
+ }
+ v = v + (it->second).shift;
+ return ;
+ }
+ }
+
dofTraits<T>::setToZero(v);
}
+
inline dataVec getDofValue(dataVec &value, MVertex *v, int iComp, int iField) const
{
getDofValue(value, v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
@@ -192,7 +233,6 @@ class dofManager{
inline void assemble(const Dof &R, const Dof &C, const dataMat &value)
{
if (!_current->isAllocated()) _current->allocate(unknown.size());
-
std::map<Dof, int>::iterator itR = unknown.find(R);
if (itR != unknown.end()){
std::map<Dof, int>::iterator itC = unknown.find(C);
@@ -207,10 +247,16 @@ class dofManager{
dofTraits<T>::mult(tmp , value, itFixed->second);
dofTraits<T>::neg(tmp);
_current->addToRightHandSide(itR->second, tmp);
- }
+ }else assembleLinConst(R, C, value);
}
}
+ if (itR == unknown.end())
+ {
+ assembleLinConst(R, C, value);
+ }
}
+
+
inline void assemble(std::vector<Dof> &R, std::vector<Dof> &C, fullMatrix<dataMat> &m)
{
if (!_current->isAllocated()) _current->allocate(unknown.size());
@@ -241,10 +287,17 @@ class dofManager{
dofTraits<T>::mult(tmp , m(i,j), itFixed->second);
dofTraits<T>::neg(tmp);
_current->addToRightHandSide(NR[i], tmp);
- }
+ }else assembleLinConst(R[i], C[j], m(i,j));
}
}
}
+ else
+ {
+ for (unsigned int j = 0; j < C.size(); j++)
+ {
+ assembleLinConst(R[i], C[j], m(i, j));
+ }
+ }
}
}
@@ -264,6 +317,19 @@ class dofManager{
{
_current->addToRightHandSide(NR[i], m(i));
}
+ else
+ {
+ typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ itConstraint = constraints.find(R[i]);
+ if (itConstraint != constraints.end())
+ {
+ for (int i=0;i<(itConstraint->second).linear.size();i++)
+ {
+ std::map<Dof, int>::iterator itC = unknown.find((itConstraint->second).linear[i].first); // lin dep in unknown ?!
+ _current->addToRightHandSide(itC->second, m(i)*(itConstraint->second).linear[i].second);
+ }
+ }
+ }
}
}
@@ -297,10 +363,16 @@ class dofManager{
dofTraits<T>::mult(tmp , m(i,j), itFixed->second);
dofTraits<T>::neg(tmp);
_current->addToRightHandSide(NR[i], tmp);
- }
+ } else assembleLinConst(R[i],R[j],m(i,j));
}
}
}
+ else
+ {
+ for (unsigned int j = 0; j < R.size(); j++){
+ assembleLinConst(R[i],R[j],m(i,j));
+ }
+ }
}
}
@@ -359,5 +431,41 @@ class dofManager{
else
return 0;
}
+
+ inline void setLinearConstraint (Dof key, DofAffineConstraint<dataVec> &affineconstraint)
+ {
+ constraints.insert(std::make_pair(key,affineconstraint));
+ }
+
+ inline void assembleLinConst(const Dof &R, const Dof &C, const dataMat &value)
+ {
+ std::map<Dof, int>::iterator itR = unknown.find(R);
+ if (itR != unknown.end())
+ {
+ typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ itConstraint = constraints.find(C);
+ if (itConstraint != constraints.end())
+ {
+ for (int i=0;i<(itConstraint->second).linear.size();i++)
+ {
+ assemble(R,(itConstraint->second).linear[i].first,(itConstraint->second).linear[i].second*value);
+ }
+ _current->addToRightHandSide(itR->second, -value * (itConstraint->second).shift);
+ }
+ }
+ else // test function ; (no shift ?)
+ {
+ typename std::map<Dof,DofAffineConstraint<dataVec> >::iterator itConstraint;
+ itConstraint = constraints.find(R);
+ if (itConstraint != constraints.end())
+ {
+ for (int i=0;i<(itConstraint->second).linear.size();i++)
+ {
+ assemble((itConstraint->second).linear[i].first,C,(itConstraint->second).linear[i].second*value);
+ }
+ }
+ }
+ }
+
};
#endif
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 91fa27586d665bcb57ce30e19e2b9db84d70084c..ac42f08dc65057293e7d3b79957b2d9c52e54022 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -370,7 +370,7 @@ class xFemFunctionSpace : public FunctionSpace<T>
FunctionSpace<T>* _spacebase;
//Function<double>* _enrichment;
public:
- virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess);
+ virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess){};
xFemFunctionSpace(FunctionSpace<T>* spacebase) : _spacebase(spacebase) {};
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){};
@@ -396,12 +396,12 @@ class FilteredFunctionSpace : public FunctionSpace<T>
public:
- virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess);
+ virtual void hessfuvw(MElement *ele, double u, double v, double w,std::vector<HessType> &hess){};
FilteredFunctionSpace<T,F>(FunctionSpace<T>* spacebase,F * filter) : _spacebase(spacebase),_filter(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 int getNumKeys(MElement *ele){};
virtual void getKeys(MElement *ele, std::vector<Dof> &keys){};
};
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index 19f21393c616a547883b23d82f71063e1507c8f0..a7f3f43e68de0b23d384b1dff58fa40bea5506f7 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -180,7 +180,37 @@ template<class Iterator,class Assembler> void NumberDofs(FunctionSpaceBase &spac
space.getKeys(e,R);
int nbdofs=R.size();
for (int i=0;i<nbdofs;++i) assembler.numberDof(R[i]);
+
}
}
+//// Mean HangingNodes
+//template <class Assembler> void FillHangingNodes(FunctionSpaceBase &space,std::map<int,std::vector <int> > &HangingNodes, Assembler &assembler, int &field, int &dim)
+//{
+// std::map<int,std::vector <int> >::iterator ith;
+// ith = HangingNodes.begin();
+// int compt = 1;
+// while (ith!=HangingNodes.end())
+// {
+// float fac;
+// fac = 1.0/(ith->second).size();
+// for (int j=0;j<dim;j++)
+// {
+// DofAffineConstraint<double> constraint;
+// int type = Dof::createTypeWithTwoInts(j, field);
+// Dof hgnd(ith->first,type);
+// for (int i=0;i<(ith->second).size();i++)
+// {
+// Dof key((ith->second)[i],type);
+// std::pair<Dof, double > linDof(key,fac);
+// constraint.linear.push_back(linDof);
+// }
+// constraint.shift = 0;
+// assembler.setLinearConstraint (hgnd, constraint);
+// }
+// ith++;
+// compt++;
+// }
+//}
+
#endif// _SOLVERALGORITHMS_H_
diff --git a/Solver/terms.h b/Solver/terms.h
index 893c9d9cef795999dd05b7721db2fcd859ed5509..b110de521a9cbedc8d727f1834e50513b5eea5bd 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -214,10 +214,10 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
fullMatrix<double> B(6, nbFF);
fullMatrix<double> BTH(nbFF, 6);
fullMatrix<double> BT(nbFF, 6);
- printf("npts : %d\n",npts);
+ //printf("npts : %d\n",npts);
m.resize(nbFF, nbFF);
m.setAll(0.);
- std::cout << m.size1() << " " << m.size2() << std::endl;
+ //std::cout << m.size1() << " " << m.size2() << std::endl;
for (int i = 0; i < npts; i++)
{
const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
@@ -235,6 +235,7 @@ class IsotropicElasticTerm : public BilinearTerm<SVector3,SVector3>
BT(j, 5) = B(5, j) = Grads[j](0,2)+Grads[j](2,0);
}
BTH.setAll(0.);
+
BTH.gemm(BT, H);
m.gemm(BTH, B, weight * detJ, 1.);
}
diff --git a/contrib/arc/CMakeLists.txt b/contrib/arc/CMakeLists.txt
deleted file mode 100644
index 2279812a9ebed9e31cd6ea4270f1db3ccae5ea9a..0000000000000000000000000000000000000000
--- a/contrib/arc/CMakeLists.txt
+++ /dev/null
@@ -1,8 +0,0 @@
-# pour boris:
-## boris cela ne marche que si le fichier est dans SVN ...
-#include(/home/boris/MyGmshProjects/CMakeLists.txt)
-# tests pour eric et christophe:
-list(APPEND EXTERNAL_INCLUDES contrib/arc)
-add_executable(elastic EXCLUDE_FROM_ALL contrib/arc/mainElasticity.cpp ${GMSH_SRC})
-target_link_libraries(elastic ${LINK_LIBRARIES})
-
diff --git a/contrib/arc/Classes/FuncGradDisc.h b/contrib/arc/Classes/FuncGradDisc.h
deleted file mode 100644
index dc09e62df6ac4c9e32351433c58b7319968337b7..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/FuncGradDisc.h
+++ /dev/null
@@ -1,276 +0,0 @@
-//
-// 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 {
-
-
- // --- F2 --- //
-
- 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++)
- {
- //std::cout<<"val[i] :" << val[i] << "\n";
- //std::cout<<"ls(i) :" << (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()) << "\n";
- f = f + std::abs((*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z()))*val[i];
- }
- f = f - std::abs((*_ls)(x,y,z));
-
- //std::cout<<"val f :" << f << "\n";
- return f;
-
-
- // --- F1 --- //
-
-
-// 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 + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*val[i];
-// }
-// f = std::abs(f);
-// 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
- {
-
-
- // ---- F2 ---- //
-
- 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 + std::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 + std::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 + std::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 + std::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 + std::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 + std::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;
- }
- }
- }
-
-
- // ---- F1 ------ //
-
-//
-// 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 + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdx;
-// dfdy = dfdy + (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdy;
-// 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 - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdx;
-// dfdy = dfdy - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdy;
-// dfdz = dfdz - (*_ls)(e->getVertex(i)->x(),e->getVertex(i)->y(),e->getVertex(i)->z())*dNdz;
-// }
-// }
-// }
-
-
-};
-
-#endif
diff --git a/contrib/arc/Classes/FuncHeaviside.h b/contrib/arc/Classes/FuncHeaviside.h
deleted file mode 100644
index 4cd9895d362bd42f074482ef61660a52de2dfdfc..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/FuncHeaviside.h
+++ /dev/null
@@ -1,39 +0,0 @@
-//
-// 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/ImageSolver.cpp b/contrib/arc/Classes/ImageSolver.cpp
deleted file mode 100644
index 8bb0a810198c4b3b6403904f4f3aaf87413876ec..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/ImageSolver.cpp
+++ /dev/null
@@ -1,377 +0,0 @@
-//
-// 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 "ImageSolver.h"
-#include "DILevelset.h"
-#include "NodeEnrichedFS.cpp"
-#include "gLevelSetMesh.h"
-
-#include "groupOfElements.h"
-#include "FuncGradDisc.h"
-
-#include "xFemFunctionSpace.cpp"
-
-
-ImageSolver::~ImageSolver()
-{
- //delete _funcEnrichment;
-}
-
-
-void ImageSolver::setMesh(GModel *m, int dim, std::vector< float > &ListLSValue, std::map< int ,std::vector<int> > & ListHangingNodes)
-{
- _dim = dim;
- int tag = 1;
- _tag = tag;
- //if (_ls) delete _ls;
- std::vector<float>::iterator itl;
- itl = ListLSValue.begin();
- int k = 1;
- while (itl!=ListLSValue.end()){
- _LevelSetValue[k] = ListLSValue[k-1] ;
- k++;
- itl++;
- }
- gLevelSetMesh *ls = new gLevelSetMesh(&_LevelSetValue,m,tag);
- GModel *pcut = m->buildCutGModel(ls);
- pModel = pcut;
- std::string ModelName = "Tree_cutLS.msh" ;
- pcut->writeMSH(ModelName,2.1,false,false);
- GModel::setCurrent(pModel);
- _ls = new gLevelSetMesh(&_LevelSetValue,m,tag);
- //ListHangingNodes.clear();
- _ListHangingNodes = & ListHangingNodes;
- //delete p;
- delete ls;
-}
-
-
-void ImageSolver::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 {
- Msg::Error("Invalid input : %s", what);
-// return;
- }
- }
- fclose(f);
-}
-
-
-void ImageSolver::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);
- }
-
- // fill mean hanging nodes
-
- std::map<int,std::vector <int> >::iterator ith;
- ith = _ListHangingNodes->begin();
- int compt = 1;
- while (ith!=_ListHangingNodes->end())
- {
- float fac;
- fac = 1.0/(ith->second).size(); // mean hanging nodes
- for (int j=0;j<_dim;j++)
- {
- DofAffineConstraint<double> constraint;
- int type = Dof::createTypeWithTwoInts(j, _tag);
- Dof hgnd(ith->first,type);
- for (int i=0;i<(ith->second).size();i++)
- {
- Dof key((ith->second)[i],type);
- std::pair<Dof, double > linDof(key,fac);
- constraint.linear.push_back(linDof);
- }
- constraint.shift = 0;
- pAssembler->setLinearConstraint (hgnd, constraint);
- }
- ith++;
- compt++;
- }
-
-// std::map<int , std::vector < int > >::iterator ith;
-// ith = _ListHangingNodes->begin();
-// std::cout << "HangingNodes :"<< std::endl;
-// while (ith!=_ListHangingNodes->end())
-// {
-//
-// std::cout<< ith->first << " " << (ith->second)[0] << " " << (ith->second)[1] <<"\n";
-// MVertex *v1;
-// MVertex *v2;
-// MVertex *v3;
-// v1 = pModel->getMeshVertexByTag(ith->first);
-// v2 = pModel->getMeshVertexByTag((ith->second)[0]);
-// v3 = pModel->getMeshVertexByTag((ith->second)[1]);
-// std::cout<<"xyz : "<<v1->x()<<" "<< v1->y() <<" "<< v1->z() << " :: " <<v2->x()<<" "<< v2->y() <<" "<< v2->z() << " :: " <<v3->x()<<" "<< v3->y() <<" "<< v3->z() << " \n ";
-// ith++;
-// }
-
-
-// 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);
-
-}
-
-
-
-PView* ImageSolver::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/ImageSolver.h b/contrib/arc/Classes/ImageSolver.h
deleted file mode 100644
index c1de3ab7eea56d0f3d40e70f6e06492fd683f8e5..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/ImageSolver.h
+++ /dev/null
@@ -1,54 +0,0 @@
-//
-// Description :
-//
-//
-// Author: <Boris Sedji>, 01/2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-
-#ifndef _IMAGE_SOLVER_H_
-#define _IMAGE_SOLVER_H_
-
-#include "elasticitySolver.h"
-#include "simpleFunction.h"
-//#include "QuadtreeLSImage.h"
-
-#include "gLevelSetMesh.cpp"
-#include "FuncGradDisc.h"
-
-class ImageSolver : 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;
- // hanging nodes
- std::map< int ,std::vector<int> > *_ListHangingNodes;
- // level set
- gLevelSetMesh *_ls;
-
- public :
-
- ImageSolver(int tag) : elasticitySolver(tag) {}
- ~ImageSolver();
- //
- 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);
- void setMesh(GModel *m, int dim, std::vector< float > &ListLSValue, std::map< int ,std::vector<int> > & ListHangingNodes);
- // 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/NodeEnrichedFS.cpp b/contrib/arc/Classes/NodeEnrichedFS.cpp
deleted file mode 100644
index 7b82abed8725970be3cb4ec396ef85ae32a6b653..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/NodeEnrichedFS.cpp
+++ /dev/null
@@ -1,195 +0,0 @@
-//
-// Description : From function space to node enriched composant function space
-//
-//
-// Author: <Eric Bechet>::<Boris Sedji>, 01/2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-#include "FuncGradDisc.h"
-#include "NodeEnrichedFS.h"
-
-
-
-template <class T> void NodeEnrichedFS<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;
- _SupportFS->f(elep,u,v,w,valsd);
- int normaldofs=valsd.size();
-
- // nbdofs element calcul
- int nbdofs = normaldofs;
- std::vector<int> EnrichedVertex;
- 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
- }
- }
-
- int curpos=vals.size();
- vals.reserve(curpos+nbdofs);
-
- // first normal dofs
- for (int i=0;i<normaldofs;i++) vals.push_back(valsd[i]);
-
- // then enriched dofs so the order is ex:(u1x,u2x,u3x,u1y,u2y,u3y,a2x,a2y,a3x,a3y)
- if (nbdofs>normaldofs) // if enriched
- {
- // Enrichment function calcul
- SPoint3 p;
- elep->pnt(u, v, w, p);
- double func;
- func = (*_funcEnrichment)(p.x(), p.y(), p.z(),elep);
-
- for (int i=0 ;i<EnrichedVertex.size();i++)
- {
- for (int j=0;j<_EnrichComp->size();j++) vals.push_back(valsd[EnrichedVertex[i]+(_EnrichComp->at(j))*elep->getNumVertices()]*func);
- }
- }
-
-
-}
-
-template <class T> void NodeEnrichedFS<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;
-
- std::vector<GradType> gradsd;
- _SupportFS->gradf(elep,u,v,w,gradsd);
-
- int normaldofs=gradsd.size();
- int nbdofs = normaldofs;
-
- GradType val;
-
- std::vector<int> EnrichedVertex;
-
- 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
- }
- }
-
-
- std::vector<ValType> valsd;
- _SupportFS->f(elep,u,v,w,valsd);
-
- int curpos=grads.size();
- grads.reserve(curpos+nbdofs);
-
- // first normal dofs
- for (int i=0;i<normaldofs;i++) grads.push_back(gradsd[i]);
-
- // then enriched dofs so the order is ex:(u1x,u2x,u3x,u1y,u2y,u3y,a2x,a2y,a3x,a3y)
- if (nbdofs>normaldofs) // 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 gradfunc(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<EnrichedVertex.size();i++)
- {
- for (int j=0;j<_EnrichComp->size();j++)
- {
- GradType GradFunc;
- tensprod(valsd[EnrichedVertex[i]+(_EnrichComp->at(j))*elep->getNumVertices()],gradfunc,GradFunc);
- grads.push_back(gradsd[EnrichedVertex[i]+(_EnrichComp->at(j))*elep->getNumVertices()]*func + GradFunc);
- }
- }
- }
-}
-
-
-template <class T> int NodeEnrichedFS<T>::getNumKeys(MElement *ele)
-{
- MElement * elep;
- if (ele->getParent()) elep = ele->getParent();
- else elep = ele;
-
- int normaldofs = _SupportFS->getNumKeys(ele);
- int nbdofs = normaldofs;
-
- for (int i=0 ;i<elep->getNumVertices();i++)
- {
- std::set<int>::iterator it;
- it = _TagEnrichedVertex->find(elep->getVertex(i)->getNum());
- if (it!=_TagEnrichedVertex->end())
- {
- nbdofs = nbdofs + 1*_EnrichComp->size(); // enriched dof
- }
- }
- return nbdofs;
-}
-
-template <class T> void NodeEnrichedFS<T>::getKeys(MElement *ele, std::vector<Dof> &keys)
-{
- MElement * elep;
- if (ele->getParent()) elep = ele->getParent();
- else elep = ele;
-
- int normalk=_SupportFS->getNumKeys(elep);
-
- std::vector<Dof> bufk;
- bufk.reserve(normalk);
- _SupportFS->getKeys(elep,bufk);
- int normaldofs=bufk.size();
- int nbdofs = normaldofs;
-
- int curpos=keys.size();
-
- std::vector<int> EnrichedVertex;
-
- 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
- }
- }
-
- keys.reserve(curpos+nbdofs);
-
- // get keys so the order is ex:(u1x,u2x,u3x,u1y,u2y,u3y,a2x,a2y,a3x,a3y)
-
- for (int i=0;i<bufk.size();i++) keys.push_back(bufk[i]);
-
- for (int i=0 ;i<EnrichedVertex.size();i++)
- {
- for (int j=0;j<_EnrichComp->size();j++)
- {
- int i1,i2;
- Dof::getTwoIntsFromType(bufk[EnrichedVertex[i]+_EnrichComp->at(j)*elep->getNumVertices()].getType(), i1,i2);
- keys.push_back(Dof(bufk[EnrichedVertex[i]+_EnrichComp->at(j)*elep->getNumVertices()].getEntity(),Dof::createTypeWithTwoInts(_EnrichComp->at(j),i1+1)));
- }
- }
-}
diff --git a/contrib/arc/Classes/NodeEnrichedFS.h b/contrib/arc/Classes/NodeEnrichedFS.h
deleted file mode 100644
index b32ca902e719cdedc7a3c307ff790cf56a86590c..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/NodeEnrichedFS.h
+++ /dev/null
@@ -1,69 +0,0 @@
-//
-// Description : From function space to node enriched composant function space
-//
-//
-// Author: <Eric Bechet>::<Boris Sedji>, 01/2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-#ifndef _NODEENRICHEDFS_H_
-#define _NODEENRICHEDFS_H_
-
-#include "SPoint3.h"
-#include "SVector3.h"
-#include "STensor3.h"
-#include <vector>
-#include <iterator>
-#include <iostream>
-#include "Numeric.h"
-#include "MElement.h"
-#include "MVertex.h"
-#include "MPoint.h"
-#include "dofManager.h"
-#include "functionSpace.h"
-#include "simpleFunction.h"
-#include "FuncGradDisc.h"
-#include <set>
-
-
-template<class T>
-class NodeEnrichedFS : public FunctionSpace<T>
-{
-
-public:
-
- typedef typename TensorialTraits<T>::ValType ValType;
- typedef typename TensorialTraits<T>::GradType GradType;
-
-
-protected:
-
- std::set<int > *_TagEnrichedVertex;
- std::vector<int> * _EnrichComp;
- //simpleFunction<double> *_funcEnrichment;
- FuncGradDisc *_funcEnrichment;
- FunctionSpace<T> *_SupportFS;
-
-public:
-
- NodeEnrichedFS(FunctionSpace<T> * SupportFS, std::set<int > * TagEnrichedVertex , std::vector<int> * EnrichComp, FuncGradDisc *funcEnrichment)
- {
- _SupportFS = SupportFS;
- _TagEnrichedVertex = TagEnrichedVertex;
- _funcEnrichment = funcEnrichment;
- _EnrichComp = EnrichComp;
- }
- virtual ~NodeEnrichedFS() {}
- // Shape function value in element at uvw
- virtual void f(MElement *ele, double u, double v, double w, 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
- virtual int getNumKeys(MElement *ele);
- // Dof keys for the element
- virtual void getKeys(MElement *ele, std::vector<Dof> &keys) ;
-};
-
-#endif
diff --git a/contrib/arc/Classes/OctreeLSImage.cpp b/contrib/arc/Classes/OctreeLSImage.cpp
deleted file mode 100644
index e248d77187288686c7ac2bfd54f83a89a036ad5a..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/OctreeLSImage.cpp
+++ /dev/null
@@ -1,1490 +0,0 @@
-//
-// Description:
-//
-//
-// Author: <Boris Sedji>, 12/2009
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-#include "GModel.h"
-#include "OctreeLSImage.h"
-#include "MElement.h"
-
-OctreeLSImage::OctreeLSImage(itk::Image< float, 3 >::Pointer image){
-
- itk::Image< float, 3 >::RegionType region;
- // Define a region to get size of the image
- region = image->GetLargestPossibleRegion ();
- // Size will be a power of two larger than the image size
- int size[3];
- size[0] = 1;
- size[1] = 1;
- size[2] = 1;
-
- _ImageSize[0] = region.GetSize(0);
- _ImageSize[1] = region.GetSize(1);
- _ImageSize[2] = region.GetSize(2);
-
- while (size[0]<_ImageSize[0]){
- size[0] = size[0]*2;
- }
- while (size[1]<_ImageSize[1]){
- size[1] = size[1]*2;
- }
- while (size[2]<_ImageSize[2]){
- size[2] = size[2]*2;
- }
-
- BoxData *data = new BoxData;
-
- if (size[0]>= size[1] & size[0]>=size[2] )
- {
- // The window of the image root
- data->boundingbox[0]=0;
- data->boundingbox[1]=size[0];
- data->boundingbox[2]=0;
- data->boundingbox[3]=size[0];
- data->boundingbox[4]=0;
- data->boundingbox[5]=size[0];
- }
- if (size[1]>= size[0] & size[1]>=size[2] )
- {
- // The window of the image root
- data->boundingbox[0]=0;
- data->boundingbox[1]=size[1];
- data->boundingbox[2]=0;
- data->boundingbox[3]=size[1];
- data->boundingbox[4]=0;
- data->boundingbox[5]=size[1];
- }
-
- if (size[2]>= size[1] & size[2]>=size[0] )
- {
- // The window of the image root
- data->boundingbox[0]=0;
- data->boundingbox[1]=size[2];
- data->boundingbox[2]=0;
- data->boundingbox[3]=size[2];
- data->boundingbox[4]=0;
- data->boundingbox[5]=size[2];
- }
-
- // initialisation of the octree
- _root = new Box(data,NULL,this);
- _root->FillLevelSetValue(image,data);
- _image = image;
- _LeafNumber = 0;
-
- std::cout<<"size root :"<<_root->BoxSize()<<"\n";
-
-}
-
-
-void OctreeLSImage::Mesh(int maxsize,int minsize){
-
-
- if (this->_root->HaveChildren()){
- return;
- }
- if (_root->BoxSize()>maxsize || !_root->IsHomogeneous()){
- _root->Mesh( maxsize, minsize);
- };
-//
-// if(_root->BoxSize()>=maxsize){
-// _root->Mesh( maxsize, minsize);
-// };
-
-}
-
-
-Box::Box(BoxData* data, Box* parent,OctreeLSImage* octree){
-
- _data = data;
- _parent = parent;
- _children[0] = NULL;
- _octree = octree;
-
-}
-
-bool Box::HaveChildren(){
-
- if (this->_children[0]==NULL) return false;
- else return true;
-
-}
-
-
-void Box::FillLevelSetValue(itk::Image< float, 3 >::Pointer image, BoxData *data){
-
- float pixelValue;
- itk::Image< float, 3 >::IndexType pixelIndex;
- itk::Image< float, 3 >::RegionType region;
- region = image->GetLargestPossibleRegion ();
-
- // fill with pixel value, and if out of region size, fill with value at largest region size
- for (int i=0;i<2;i++){
- for (int j=0;j<2;j++){
- for (int k=0;k<2;k++){
-
- if (data->boundingbox[i] < region.GetSize(0)) pixelIndex[0] = data->boundingbox[i];
- else pixelIndex[0] = region.GetSize(0)-1; // x position
-
- if (data->boundingbox[2+j] < region.GetSize(1)) pixelIndex[1] = data->boundingbox[2+j];
- else pixelIndex[1] = region.GetSize(1)-1; // y position
-
- if (data->boundingbox[4+k] < region.GetSize(2)) pixelIndex[2] = data->boundingbox[4+k];
- else pixelIndex[2] = region.GetSize(2)-1; // z position
-
- pixelValue = image->GetPixel( pixelIndex );
- data->LevelSetValue[i][j][k] = pixelValue;
-
- }
- }
- }
-
-}
-
-int Box::BoxSize(){
-
- int SizeX;
- int SizeY;
- int SizeZ;
- int size;
-
- SizeX = this->_data->boundingbox[1] - this->_data->boundingbox[0];
- SizeY = this->_data->boundingbox[3] - this->_data->boundingbox[2];
- SizeZ = this->_data->boundingbox[5] - this->_data->boundingbox[4];
-
- if (SizeX<SizeY) size = SizeX;
- else size = SizeY;
-
- if (SizeZ<size) size = SizeZ;
-
- return size;
-
-}
-
-bool Box::IsHomogeneous(){
-
- float OneValue = this->_data->LevelSetValue[0][0][0];
-
- for (int i=0;i<2;i++){
- for (int j=0;j<2;j++){
- for (int k=0;k<2;k++){
- OneValue = OneValue*this->_data->LevelSetValue[i][j][k];
- if (OneValue<0) {
- return false; // if sign change => not homogeneous
- };
- OneValue = this->_data->LevelSetValue[i][j][k];
- }
- }
- }
- return true;
-
-}
-
-
-// Recursive dividing function - eight nodes
-
-void Box::Mesh(int & maxsize, int & minsize){
-
- if (((this->BoxSize()<=maxsize) & (this->IsHomogeneous())) | (this->BoxSize()<=minsize)){ // (max size & homogenous) ou (min size) => dont divide
- return; // dont divide
- }
-
-// std::cout<<"databoundingbox : "<< this->_data->boundingbox[1] <<"\n";
-// if(((this->BoxSize()<=maxsize) & !(this->_data->boundingbox[0] == 0 & this->_data->boundingbox[3] == 1024 & this->_data->boundingbox[4] == 0)) | (this->BoxSize()<=minsize)){ // (max size & homogenous) ou (min size) => dont divide
-// return; // dont divide
-// }
-
- // else we divide...
- this->Divide();
-
- for (int n = 0;n<8;n++){
- _children[n]->Mesh(maxsize,minsize);
- }
-
-}
-
-void Box::PrintLevelSetValue()
-{
-
- for (int i=0;i<2;i++){
- for (int j=0;j<2;j++){
- for (int k=0;k<2;k++){
- std::cout<<"\n "<< this->_data->LevelSetValue[i][j][k];
- }
- }
- }
-
-}
-
-
-
-GModel *OctreeLSImage::CreateGModel(bool simplex, double facx, double facy, double facz,int & sizemax,int & sizemin){
-
- this->FillMeshInfo(sizemin);
- this->FillHangingNodes();
-
- std::cout<<"HangingNodes size :"<< _HangingNodes.size()<<"\n";
-
- itk::Image< float, 3 >::RegionType region;
- // Define a region to get size of the image
- region = _image->GetLargestPossibleRegion ();
-
-
- std::map<int, MVertex*> vertexMap;
- std::map<int,int> vertexNum;
- std::vector<MVertex*> vertexVector;
-
- std::vector<int> numElement;
- std::vector<std::vector<int> > vertexIndices;
- std::vector<int> elementType ;
- std::vector<int> physical;
- std::vector<int> elementary;
- std::vector<int> partition;
- std::vector<double> d;
- std::map<int, std::vector<double> > data;
- data.clear();
-
- double eps = 1e-6;
-
- int numVertices;
-
- std::vector<std::vector<int> >::iterator ite;
-
- int i = 1;
- int k = 1;
-
- int a;
-
- int xlimit = _ImageSize[0] - _ImageSize[0]%sizemax; // lost of image pixels...
- int ylimit = _ImageSize[1] - _ImageSize[1]%sizemax;
- int zlimit = _ImageSize[2] - _ImageSize[2]%sizemax;
-
- std::cout<<"\nxlimit :"<< xlimit<<"\n";
- std::cout<<"ylimit :"<< ylimit<<"\n";
- std::cout<<"zlimit :"<< zlimit<<"\n";
- std::cout<<"Listnodes size :"<<_ListNodes.size()<<"\n";
-
- i = 1;
- ite = _ListElements.begin();
- if (!simplex) // first try, to be improved
- {
- while (ite!=_ListElements.end()) {
- int num, type, phys = 0, ele = 0, part = 0, numVertices;
- num = i;
- type = 5;
- ele = 7;
- phys = 7;
- part = 0;
- numVertices = MElement::getInfoMSH(type);
- std::vector <int> indices;
- for (int j = 0; j < numVertices; j++){
- indices.push_back((*ite)[j]);
- }
- numElement.push_back(i);
- vertexIndices.push_back(indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- indices.clear();
- ite++;
- i++;
- }
- }
- else // if simplex
- {
- while (ite!=_ListElements.end()) {
-
- bool inImage;
- inImage = true;
- for (int j = 0;j<8;j++){
- if ( _ListNodes[(*ite)[j]-1][0] > xlimit | _ListNodes[(*ite)[j]-1][1] > ylimit | _ListNodes[(*ite)[j]-1][2] > zlimit ){
- inImage = false;
- }
- }
- if (inImage)
- {
- k++;
- for (int j = 0;j<8;j++)
- {
- if (vertexMap.find((*ite)[j]) == vertexMap.end()) // if not in
- {
- float xyz[3];
- d.clear();
- xyz[0] = (_ListNodes[(*ite)[j]-1][0])*facx;
- xyz[1] = (_ListNodes[(*ite)[j]-1][1])*facy;
- xyz[2] = (_ListNodes[(*ite)[j]-1][2])*facz;
- MVertex* newVertex = new MVertex(xyz[0], xyz[1], xyz[2], 0, (*ite)[j]);
- vertexMap[(*ite)[j]] = newVertex;
- d.push_back(_ListLSValue[(*ite)[j]-1]);
- data[(*ite)[j]]=d;
- d.clear();
- }
- }
-
- int num, type, phys = 0, ele = 0, part = 0, numVertices;
- num = i;
- type = 4;
- ele = 7;
- phys = 7;
- part = 0;
- numVertices = 4; //
- std::vector <int> indices;
- std::vector <int> front_1_indices;
- std::vector <int> front_2_indices;
- std::vector <int> front_3_indices;
- std::vector <int> front_4_indices;
- std::vector <int> front_5_indices;
- std::vector <int> front_6_indices;
-
- // --- Domain element tetraedre 1 ---
- indices.clear();
- indices.push_back((*ite)[0]);
- indices.push_back((*ite)[1]);
- indices.push_back((*ite)[3]);
- indices.push_back((*ite)[7]);
- type = 4;
- ele = 7;
- phys = 7;
- part = 0;
- numElement.push_back(i);
- vertexIndices.push_back(indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
-
-
- // --- frontier element---
- front_1_indices.clear();
- front_2_indices.clear();
- front_3_indices.clear();
- front_4_indices.clear();
- front_5_indices.clear();
- front_6_indices.clear();
- for (int j =0;j<4;j++)
- {
- if (vertexMap[indices[j]]->x() < eps*facx)
- {
- front_1_indices.push_back(indices[j]);
- }
- if (vertexMap[indices[j]]->y() < eps*facy)
- {
- front_2_indices.push_back(indices[j]);
- }
- }
- if (front_2_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 2;
- phys = 2;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_2_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
- if (front_1_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 1;
- phys = 1;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_1_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
- // --- Domain element tetraedre 2 ---
- indices.clear();
- indices.push_back((*ite)[1]);
- indices.push_back((*ite)[2]);
- indices.push_back((*ite)[3]);
- indices.push_back((*ite)[7]);
- type = 4;
- ele = 7;
- phys = 7;
- part = 0;
- numElement.push_back(i);
- vertexIndices.push_back(indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
-
- // --- frontier element tetraedre 2---
- front_1_indices.clear();
- front_2_indices.clear();
- front_3_indices.clear();
- front_4_indices.clear();
- front_5_indices.clear();
- front_6_indices.clear();
- for (int j =0;j<4;j++)
- {
- if (vertexMap[indices[j]]->x() < eps*facx)
- {
- front_1_indices.push_back(indices[j]);
- }
- if (vertexMap[indices[j]]->z() > (zlimit-eps)*facz)
- {
- front_6_indices.push_back(indices[j]);
- }
- }
-
- if (front_1_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 1;
- phys = 1;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_1_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
- if (front_6_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 6;
- phys = 6;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_6_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
-
- // --- Domain element tetraedre 3 ---
- indices.clear();
- indices.push_back((*ite)[1]);
- indices.push_back((*ite)[6]);
- indices.push_back((*ite)[2]);
- indices.push_back((*ite)[7]);
- type = 4;
- ele = 7;
- phys = 7;
- part = 0;
- numElement.push_back(i);
- vertexIndices.push_back(indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
-
-
- // --- frontier element tetraedre 3---
- front_1_indices.clear();
- front_2_indices.clear();
- front_3_indices.clear();
- front_4_indices.clear();
- front_5_indices.clear();
- front_6_indices.clear();
- for (int j =0;j<4;j++)
- {
- if (vertexMap[indices[j]]->y() > (ylimit-eps)*facy)
- {
- front_5_indices.push_back(indices[j]);
- }
- if (vertexMap[indices[j]]->z() > (zlimit-eps)*facz)
- {
- front_6_indices.push_back(indices[j]);
- }
- }
-
- if (front_5_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 5;
- phys = 5;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_5_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
- if (front_6_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 6;
- phys = 6;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_6_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
-
- // --- Domain element tetraedre 4 ---
- indices.clear();
- indices.push_back((*ite)[1]);
- indices.push_back((*ite)[5]);
- indices.push_back((*ite)[6]);
- indices.push_back((*ite)[7]);
- type = 4;
- type = 4;
- ele = 7;
- phys = 7;
- part = 0;
- numElement.push_back(i);
- vertexIndices.push_back(indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
-
- // --- frontier element tetraedre 4---
- front_1_indices.clear();
- front_2_indices.clear();
- front_3_indices.clear();
- front_4_indices.clear();
- front_5_indices.clear();
- front_6_indices.clear();
- for (int j =0;j<4;j++)
- {
- if (vertexMap[indices[j]]->x() > (xlimit-eps)*facx)
- {
- front_4_indices.push_back(indices[j]);
- }
- if (vertexMap[indices[j]]->y() > (ylimit-eps)*facy)
- {
- front_5_indices.push_back(indices[j]);
- }
- }
- if (front_4_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 4;
- phys = 4;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_4_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
- if (front_5_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 5;
- phys = 5;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_5_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
- // --- Domain element tetraedre 5 ---
- indices.clear();
- indices.push_back((*ite)[1]);
- indices.push_back((*ite)[5]);
- indices.push_back((*ite)[7]);
- indices.push_back((*ite)[4]);
- type = 4;
- ele = 7;
- phys = 7;
- part = 0;
- numElement.push_back(i);
- vertexIndices.push_back(indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
-
- // --- frontier element tetraedre 5---
- front_1_indices.clear();
- front_2_indices.clear();
- front_3_indices.clear();
- front_4_indices.clear();
- front_5_indices.clear();
- front_6_indices.clear();
- for (int j =0;j<4;j++)
- {
- if (vertexMap[indices[j]]->z() < eps*facz)
- {
- front_3_indices.push_back(indices[j]);
- }
- if (vertexMap[indices[j]]->x() > (xlimit-eps)*facx)
- {
- front_4_indices.push_back(indices[j]);
- }
- }
-
- if (front_3_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 3;
- phys = 3;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_3_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
- if (front_4_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 4;
- phys = 4;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_4_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
- // --- Domain element tetraedre 6 ---
- indices.clear();
- indices.push_back((*ite)[1]);
- indices.push_back((*ite)[4]);
- indices.push_back((*ite)[7]);
- indices.push_back((*ite)[0]);
- type = 4;
- type = 4;
- ele = 7;
- phys = 7;
- part = 0;
- numElement.push_back(i);
- vertexIndices.push_back(indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
-
-
- // --- frontier element tetraedre 6---
- front_1_indices.clear();
- front_2_indices.clear();
- front_3_indices.clear();
- front_4_indices.clear();
- front_5_indices.clear();
- front_6_indices.clear();
- for (int j =0;j<4;j++)
- {
- if (vertexMap[indices[j]]->y() < eps*facy)
- {
- front_2_indices.push_back(indices[j]);
- }
- if (vertexMap[indices[j]]->z() < eps*facz)
- {
- front_3_indices.push_back(indices[j]);
- }
- }
- if (front_2_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 2;
- phys = 2;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_2_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
- if (front_3_indices.size() == 3)
- {
- type = 2; // MSH_TRI
- ele = 3;
- phys = 3;
- part = 0;
- numVertices = 3;
- numElement.push_back(i);
- vertexIndices.push_back(front_3_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
- }
- ite++; // next octree leaf
- }
- }
-
- std::cout<<"numElement size :"<<numElement.size()<<"\n";
- std::cout<<"vertexIndices size :"<<vertexIndices.size()<<"\n";
-
- std::cout<<"data size :"<<data.size()<<"\n";
- std::cout<<"vertexMap size :"<<vertexMap.size()<<"\n";
- GModel *gmod = GModel::createGModel(vertexMap,numElement,vertexIndices,elementType,physical,elementary,partition);
- // Write a .msh file with the level set values as postview data
- std::string postTypeName = "LevelSet Value" ;
- PView *pv = new PView (postTypeName, "NodeData", gmod, data);
- //PView *pv = octree.CreateLSPView(m);
- bool useadapt = true;
- pv->getData(useadapt)->writeMSH("LSPView.msh", false);
-
- std::cout<<"getNumScalars :"<< pv->getData(useadapt)->getNumScalars()<<"\n";
-
-//
-// std::string LevelSetFileName = "LevelSetValue.msh" ;
-// FILE *f = fopen(LevelSetFileName.c_str(), "w");
-// char name[256] = "LSValues";
-// int num;
-// float val;
-// int numnodes;
-// fprintf(f, "%s\n", name);
-// fprintf(f, "%d\n", data.size());
-// i = 1;
-// while(data.find(i)!=data.end()){
-// fprintf(f,"%d %f\n",i,data[i][0]);
-// i++;
-// }
-// fclose(f);
-
- return gmod;
-
-}
-
-PView *OctreeLSImage::CreateLSPView(GModel* m){
-
- std::map<int, std::vector<double> > data;
- std::vector<float>::iterator itf;
- itf = _ListLSValue.begin();
- int i = 1;
-
- while (itf!=_ListLSValue.end()){
- std::vector<double> d;
- d.push_back((*itf)) ;
- data[i] = d;
- d.clear();
- itf++;
- i++;
- }
- std::string postTypeName = "LevelSet Value" ;
- PView *pv = new PView (postTypeName, "NodeData", m, data, 0.0);
-
- return pv;
-
-}
-
-
-void OctreeLSImage::FillHangingNodes(){
-
- _root->FillHangingNodes(_ListNodes,_ListElements,_HangingNodes);
-
-}
-
-
-void Box::FillHangingNodes(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements,std::map< int ,std::vector<int> > &HangingNodes){
-
- std::vector<std::vector<int> >::iterator it;
-
- it = ListNodes.begin();
- int k=0;
- while (it!=ListNodes.end())
- {
- std::vector<Box*> Leafs;
- std::vector<int> MasterNodes;
- Leafs.clear();
- MasterNodes.clear();
- GetLeafsWith(Leafs, (*it)[0], (*it)[1],(*it)[2]);
- if (Leafs.size()==5)
- {
- Box *MaxLeaf;
- MaxLeaf = Leafs[0];
- int Face[4];
- for (int i=1;i<5;i++)
- {
- if (MaxLeaf->BoxSize() < Leafs[i]->BoxSize())
- MaxLeaf = Leafs[i];
- }
- int j;
- j = 0;
- for (int i = 0 ; i<8;i++)
- {
- if (ListNodes[MaxLeaf->_ElementNodes[i]-1][0] == (*it)[0] | ListNodes[MaxLeaf->_ElementNodes[i]-1][1]== (*it)[1] | ListNodes[MaxLeaf->_ElementNodes[i]-1][2]== (*it)[2])
- {
- // std::cout<<"NodeIn xyz : "<< ListNodes[MaxLeaf->_ElementNodes[i]-1][0] << " " << ListNodes[MaxLeaf->_ElementNodes[i]-1][1] << " " << ListNodes[MaxLeaf->_ElementNodes[i]-1][2] <<"\n";
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- if (MasterNodes.size()<5) Face[j] = i;
- j++;
- }
-
- }
- if (MasterNodes.size()==4)
- {
- MasterNodes.clear();
- if (Face[0] == 0 & Face[1] == 1 & Face[2] == 2 & Face[3] == 3 )
- {
- MasterNodes.push_back(MaxLeaf->_ElementNodes[1]);
- MasterNodes.push_back(MaxLeaf->_ElementNodes[3]);
- }
- if (Face[0] == 0 & Face[1] == 3 & Face[2] == 4 & Face[3] == 7 )
- {
- MasterNodes.push_back(MaxLeaf->_ElementNodes[0]);
- MasterNodes.push_back(MaxLeaf->_ElementNodes[7]);
- }
- if (Face[0] == 0 & Face[1] == 1 & Face[2] == 4 & Face[3] == 5 )
- {
- MasterNodes.push_back(MaxLeaf->_ElementNodes[1]);
- MasterNodes.push_back(MaxLeaf->_ElementNodes[4]);
- }
- if (Face[0] == 4 & Face[1] == 5 & Face[2] == 6 & Face[3] == 7 )
- {
- MasterNodes.push_back(MaxLeaf->_ElementNodes[5]);
- MasterNodes.push_back(MaxLeaf->_ElementNodes[7]);
- }
- if (Face[0] == 1 & Face[1] == 2 & Face[2] == 5 & Face[3] == 6 )
- {
- MasterNodes.push_back(MaxLeaf->_ElementNodes[1]);
- MasterNodes.push_back(MaxLeaf->_ElementNodes[6]);
- }
- if (Face[0] == 2 & Face[1] == 3 & Face[2] == 6 & Face[3] == 7 )
- {
- MasterNodes.push_back(MaxLeaf->_ElementNodes[2]);
- MasterNodes.push_back(MaxLeaf->_ElementNodes[7]);
- }
-
- }
- if (MasterNodes.size() > 4)
- {
- //std::cout<<"Recompt : \n";
- MasterNodes.clear();
- for (int i = 0 ; i<8;i++)
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0] == (*it)[0] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][1]== (*it)[1] | ListNodes[MaxLeaf->_ElementNodes[i]-1][2]== (*it)[2])
- {
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- //std::cout<<"mm x NodeIn xyz : "<< ListNodes[MaxLeaf->_ElementNodes[i]-1][0] << " " << ListNodes[MaxLeaf->_ElementNodes[i]-1][1] << " " << ListNodes[MaxLeaf->_ElementNodes[i]-1][2] <<"\n";
- }
- }
- else if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][1] == (*it)[1] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0]== (*it)[0] | ListNodes[MaxLeaf->_ElementNodes[i]-1][2]== (*it)[2])
- {
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- //std::cout<<"mm y NodeIn xyz : "<< ListNodes[MaxLeaf->_ElementNodes[i]-1][0] << " " << ListNodes[MaxLeaf->_ElementNodes[i]-1][1] << " " << ListNodes[MaxLeaf->_ElementNodes[i]-1][2] <<"\n";
- }
- }else if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][2] == (*it)[2] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0]== (*it)[0] | ListNodes[MaxLeaf->_ElementNodes[i]-1][1]== (*it)[1])
- {
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- // std::cout<<"mm z NodeIn xyz : "<< ListNodes[MaxLeaf->_ElementNodes[i]-1][0] << " " << ListNodes[MaxLeaf->_ElementNodes[i]-1][1] << " " << ListNodes[MaxLeaf->_ElementNodes[i]-1][2] <<"\n";
- }
- }
- }
- }
- HangingNodes[k+1] = MasterNodes;
- }
- if (Leafs.size()==3)
- {
- Box *MaxLeaf;
- MaxLeaf = Leafs[0];
- for (int i=1;i<3;i++)
- {
- if (MaxLeaf->BoxSize() < Leafs[i]->BoxSize())
- MaxLeaf = Leafs[i];
- }
- for (int i = 0 ; i<8;i++)
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0] == (*it)[0] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][1]== (*it)[1] | ListNodes[MaxLeaf->_ElementNodes[i]-1][2]== (*it)[2])
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- } else if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][1] == (*it)[1] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0]== (*it)[0] | ListNodes[MaxLeaf->_ElementNodes[i]-1][2]== (*it)[2])
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- } else if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][2] == (*it)[2] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0]== (*it)[0] | ListNodes[MaxLeaf->_ElementNodes[i]-1][1]== (*it)[1])
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- }
- }
- HangingNodes[k+1] = MasterNodes;
- }
- if (Leafs.size()==6)
- {
- Box *MaxLeaf;
- MaxLeaf = Leafs[0];
- for (int i=1;i<6;i++)
- {
- if (MaxLeaf->BoxSize() < Leafs[i]->BoxSize())
- MaxLeaf = Leafs[i];
- }
- for (int i = 0 ; i<8;i++)
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0] == (*it)[0] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][1]== (*it)[1] | ListNodes[MaxLeaf->_ElementNodes[i]-1][2]== (*it)[2])
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- } else if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][1] == (*it)[1] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0]== (*it)[0] | ListNodes[MaxLeaf->_ElementNodes[i]-1][2]== (*it)[2])
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- } else if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][2] == (*it)[2] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0]== (*it)[0] | ListNodes[MaxLeaf->_ElementNodes[i]-1][1]== (*it)[1])
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- }
- }
- HangingNodes[k+1] = MasterNodes;
- }
-
- if (Leafs.size()==7)
- {
- Box *MaxLeaf;
- MaxLeaf = Leafs[0];
- for (int i=1;i<7;i++)
- {
- if (MaxLeaf->BoxSize() < Leafs[i]->BoxSize())
- MaxLeaf = Leafs[i];
- }
- for (int i = 0 ; i<8;i++)
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0] == (*it)[0] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][1]== (*it)[1] | ListNodes[MaxLeaf->_ElementNodes[i]-1][2]== (*it)[2])
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- } else if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][1] == (*it)[1] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0]== (*it)[0] | ListNodes[MaxLeaf->_ElementNodes[i]-1][2]== (*it)[2])
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- } else if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][2] == (*it)[2] )
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0]== (*it)[0] | ListNodes[MaxLeaf->_ElementNodes[i]-1][1]== (*it)[1])
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- }
- }
- HangingNodes[k+1] = MasterNodes;
- }
- it++;
- k++;
- }
-}
-
-void OctreeLSImage::FillMeshInfo(int & pas){
-
- if (!_ListNodes.empty()){
- return;
- }
-
- _root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue,pas);
-
-}
-
-
-
-// fonction par balayage de région
-
-
-void Box::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector<float> &ListLSValue,int &pas){
-
-
- int size = _octree->GetRoot()->BoxSize();
-
- int xpas =pas; // Ã modifier si sizeZ n'est pas le plus petit
- int ypas = pas;
- int zpas = pas;
-
- std::vector<std::vector<int> >::iterator it;
-
- it = ListNodes.begin();
- int iti=0;
-
- for (int z = 0;z<=size;z+=zpas)
- for (int y = 0;y<=size;y+=ypas)
- for (int x = 0;x<=size; x+=xpas)
- {
- std::vector<Box*> Leafs;
- Leafs.clear();
- GetLeafsWith(Leafs, x, y, z);
- std::map<int,std::vector< int > > ElementsNodes;
- std::vector<int> NodesIn;
-
- bool added = false;
- for (unsigned int l=0;l<Leafs.size();l++){
- if (!added){
- std::vector<int> XYZ;
- XYZ.push_back(x);
- XYZ.push_back(y);
- XYZ.push_back(z);
- for (int i = 0 ; i<2;i++ )
- for (int j = 0 ; j<2;j++ )
- for (int k = 0 ; k<2;k++ ){
- if (x == Leafs[l]->GetData()->boundingbox[i] & y == Leafs[l]->GetData()->boundingbox[j+2] & z == Leafs[l]->GetData()->boundingbox[4+k] ){
- ListLSValue.push_back(Leafs[l]->GetData()->LevelSetValue[i][j][k]);
- ListNodes.push_back(XYZ);
- added = true;
- iti++;
- }
- }
- }else break;
- }
- if (added){
- for (unsigned int l=0;l<Leafs.size();l++){
- for (int i = 0 ; i<2;i++ )
- for (int k = 0 ; k<2;k++ )
- for (int j = 0 ; j<2;j++ ){
- int pos ;
- pos = i*4+k*2+j*1; // pour l'autre méthode
- if (x == Leafs[l]->GetData()->boundingbox[i] & y == Leafs[l]->GetData()->boundingbox[j+2] & z == Leafs[l]->GetData()->boundingbox[k+4]){
- Leafs[l]->SetElementNode(pos,iti); //différence avec l'autre méthode
- }
- }
- }
- }
- }
- std::cout<<"\nNumber Nodes : "<< ListNodes.size() << "\n";
- _octree->GetRoot()->FillElementsNode(ListElements);
-
-}
-
-void Box::FillElementsNode(std::vector< std::vector<int> > &ListElements){
- // if it s a leaf then fill list node if node dont exist in list
- if (!this->HaveChildren()) {
- std::vector<int> ElementNodes;
- for (int i = 0 ; i < 8 ; i++) ElementNodes.push_back(_ElementNodes[i]);
- int temp;
- temp = ElementNodes[2];
- ElementNodes[2] = ElementNodes[3];
- ElementNodes[3] = temp;
- temp = ElementNodes[6];
- ElementNodes[6] = ElementNodes[7];
- ElementNodes[7] = temp;
-
-
- temp = _ElementNodes[2];
- _ElementNodes[2] = _ElementNodes[3];
- _ElementNodes[3] = temp;
- temp = _ElementNodes[6];
- _ElementNodes[6] = _ElementNodes[7];
- _ElementNodes[7] = temp;
-
-// temp = ElementNodes[1];
-// ElementNodes[1] = ElementNodes[2];
-// ElementNodes[2] = temp;
-// temp = ElementNodes[2];
-// ElementNodes[2] = ElementNodes[3];
-// ElementNodes[3] = temp;
-//
-// temp = ElementNodes[5];
-// ElementNodes[5] = ElementNodes[6];
-// ElementNodes[6] = temp;
-// temp = ElementNodes[6];
-// ElementNodes[6] = ElementNodes[7];
-// ElementNodes[7] = temp;
-
- ListElements.push_back(ElementNodes);
-// std::cout<<"\n Nodes :";
-// for (int i = 0 ; i < 8 ; i++)
-// std::cout<<" " << ElementNodes[i];
- return;
- }
-
- for (int n = 0 ; n < 8 ; n++){
- this->_children[n]->FillElementsNode(ListElements);
- }
-}
-
-
-/*
-void Box::FillInfo(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector<float> &ListLSValue){
-
- // if it s a leaf then fill list node if node dont exist in list
- if (!this->HaveChildren()) {
- std::vector<int> NodesIn;
- NodesIn.clear();
- int temp;
- for (int i = 0 ; i<2;i++ ){
- for (int k = 0 ; k<2;k++ )
- for (int j = 0 ; j<2;j++ ){
- std::vector<int> XYZ;
- float LSValue;
- XYZ.push_back(this->GetData()->boundingbox[i]);
- XYZ.push_back(this->GetData()->boundingbox[2+j]);
- XYZ.push_back(this->GetData()->boundingbox[4+k]);
- LSValue = this->GetData()->LevelSetValue[i][j][k];
- std::vector<std::vector<int> >::iterator it;
- it = ListNodes.begin();
- bool add = true;
- int l;
- l = 1;
- while (it!=ListNodes.end()){ // attention segmentation & pas opti
- if (((*it)[0] == XYZ[0]) & ((*it)[1] == XYZ[1]) & ((*it)[2] == XYZ[2])){
- add = false;
- NodesIn.push_back(l);
- break;
- }
- it++;
- l++;
- }
- if (add) {
- ListNodes.push_back(XYZ);
- ListLSValue.push_back(LSValue);
- NodesIn.push_back(l);
- }
- XYZ.clear();
- }
- }
- temp = NodesIn[3];
- NodesIn[3] = NodesIn[2];
- NodesIn[2] = temp;
- temp = NodesIn[7];
- NodesIn[7] = NodesIn[6];
- NodesIn[6] = temp;
- ListElements.push_back(NodesIn);
- return;
- }
-
- for (int n = 0 ; n < 8 ; n++){
- this->_children[n]->FillMeshInfo(ListNodes,ListElements,ListLSValue);
- }
-
-}*/
-
-/*
-std::vector< std::vector <int> > *OctreeLSImage::GetListElements(){
-
- if (!_ListElements.empty()){
- return &_ListElements;
- }
-
- _root->GetLeafElements(_ListElements);
-
- return &_ListElements;
-
-}*/
-
-
-/*
-void Box::GetLeafElements(std::vector< std::vector<int> > &ListElements){
-
-// if (!this->HaveChildren()) {
-// std::vector<std::vector<int> >::iterator it;
-// std::vector<std::vector<int> > *ListNodes;
-// ListNodes = this->_octree->GetListNodes();
-// it = ListNodes->begin();
-// std::vector<int> NodesIn;
-// int l = 1;
-// while (it!=ListNodes->end()) {
-// for (int i = 0 ; i<2;i++ ){
-// for (int j = 0 ; j<2;j++ ){
-// for (int k = 0 ; k<2;k++ ){
-// if ((_data->boundingbox[i] == (*it)[0]) & (_data->boundingbox[2+j] == (*it)[1]) & (_data->boundingbox[4+k] == (*it)[2]))
-// NodesIn.push_back(l);
-// }
-// }
-// }
-// it++;
-// l++;
-// }
-// ListElements.push_back(NodesIn);
-// return;
-// }
-//
-// for (int n = 0 ; n < 8 ; n++){
-// this->_children[n]->GetLeafElements(ListElements);
-// }
-
-}
-*/
-
-
-
-
-void OctreeLSImage::SetLeafNumber(){
- _LeafNumber = 0;
- _root->CountLeafNumber(_LeafNumber);
-}
-
-void Box::CountLeafNumber(int &LeafNumber){
-
- if (this->GetChild(0)==NULL) {
- LeafNumber = LeafNumber + 1;
- return;
- }
-
- for (int n = 0 ; n < 8 ; n++){
- this->_children[n]->CountLeafNumber(LeafNumber);
- }
-
-}
-
-// Divide function
-
-void Box::Divide(){
-
-// box divide in 8 boxes whit limite xl = (xmin + xmax)/2 ; yl = (ymin + ymax)/2 ; zl = (zmin + zmax)/2
-
- int n;
-
- // Box 1 : < xl ; < yl ; < zl
-
- n = 0;
- BoxData* data1 = new BoxData;
-
- data1->boundingbox[0]=this->_data->boundingbox[0];
- data1->boundingbox[1]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data1->boundingbox[2]=this->_data->boundingbox[2];
- data1->boundingbox[3]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
- data1->boundingbox[4]=this->_data->boundingbox[4];
- data1->boundingbox[5]=(this->_data->boundingbox[5]+this->_data->boundingbox[4])/2;
-
- this->_children[n] = new Box(data1,this,this->_octree);
- this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
-
- // Box 2 : > xl ; < yl ; < zl
-
- n = 1;
- BoxData* data2 = new BoxData;
-
- data2->boundingbox[0]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data2->boundingbox[1]=this->_data->boundingbox[1];
- data2->boundingbox[2]=this->_data->boundingbox[2];
- data2->boundingbox[3]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
- data2->boundingbox[4]=this->_data->boundingbox[4];
- data2->boundingbox[5]=(this->_data->boundingbox[5]+this->_data->boundingbox[4])/2;
-
- this->_children[n] = new Box(data2,this,this->_octree);
- this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
-
- // Box 3 : < xl ; > yl ; < zl
-
- n = 2;
- BoxData* data3 = new BoxData;
-
- data3->boundingbox[0]=this->_data->boundingbox[0];
- data3->boundingbox[1]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data3->boundingbox[2]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
- data3->boundingbox[3]=this->_data->boundingbox[3];
- data3->boundingbox[4]=this->_data->boundingbox[4];
- data3->boundingbox[5]=(this->_data->boundingbox[5]+this->_data->boundingbox[4])/2;
-
- this->_children[n] = new Box(data3,this,this->_octree);
- this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
-
- // Box 4 : > xl ; > yl ; < zl
-
- n = 3;
- BoxData* data4 = new BoxData;
-
- data4->boundingbox[0]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data4->boundingbox[1]=this->_data->boundingbox[1];
- data4->boundingbox[2]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
- data4->boundingbox[3]=this->_data->boundingbox[3];
- data4->boundingbox[4]=this->_data->boundingbox[4];
- data4->boundingbox[5]=(this->_data->boundingbox[5]+this->_data->boundingbox[4])/2;
-
- this->_children[n] = new Box(data4,this,this->_octree);
- this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
-
- // Box 5 : < xl ; < yl ; > zl
-
- n = 4;
- BoxData* data5 = new BoxData;
-
- data5->boundingbox[0]=this->_data->boundingbox[0];
- data5->boundingbox[1]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data5->boundingbox[2]=this->_data->boundingbox[2];
- data5->boundingbox[3]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
- data5->boundingbox[4]=(this->_data->boundingbox[5]+this->_data->boundingbox[4])/2;
- data5->boundingbox[5]=this->_data->boundingbox[5];
-
- this->_children[n] = new Box(data5,this,this->_octree);
- this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
-
- // Box 6 : > xl ; < yl ; > zl
-
- n = 5;
- BoxData* data6 = new BoxData;
-
- data6->boundingbox[0]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data6->boundingbox[1]=this->_data->boundingbox[1];
- data6->boundingbox[2]=this->_data->boundingbox[2];
- data6->boundingbox[3]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
- data6->boundingbox[4]=(this->_data->boundingbox[5]+this->_data->boundingbox[4])/2;
- data6->boundingbox[5]=this->_data->boundingbox[5];
-
- this->_children[n] = new Box(data6,this,this->_octree);
- this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
-
- // Box 7 : < xl ; > yl ; > zl
-
- n = 6;
- BoxData* data7 = new BoxData;
-
- data7->boundingbox[0]=this->_data->boundingbox[0];
- data7->boundingbox[1]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data7->boundingbox[2]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
- data7->boundingbox[3]=this->_data->boundingbox[3];
- data7->boundingbox[4]=(this->_data->boundingbox[5]+this->_data->boundingbox[4])/2;
- data7->boundingbox[5]=this->_data->boundingbox[5];
-
- this->_children[n] = new Box(data7,this,this->_octree);
- this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
-
- // Box 8 : > xl ; > yl ; > zl
-
- n = 7;
- BoxData* data8 = new BoxData;
-
- data8->boundingbox[0]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data8->boundingbox[1]=this->_data->boundingbox[1];
- data8->boundingbox[2]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
- data8->boundingbox[3]=this->_data->boundingbox[3];
- data8->boundingbox[4]=(this->_data->boundingbox[5]+this->_data->boundingbox[4])/2;
- data8->boundingbox[5]=this->_data->boundingbox[5];
-
- this->_children[n] = new Box(data8,this,this->_octree);
- this->_children[n]->FillLevelSetValue(_children[n]->_octree->GetImage(),_children[n]->GetData());
-}
-
-bool OctreeLSImage::Smooth(){
-
- return _root->Smooth();
-
-}
-
-
-bool Box::Smooth(){
-
-
- if (!this->HaveChildren()) {
- bool smoothed = true;
- BoxData*data = this->GetData();
- std::vector<Box*> Leafs;
- for (int i = 0 ; i<2;i++ ){
- for (int j = 0 ; j<2;j++ ){
- for (int k = 0 ; k<2;k++ ){
- Leafs.clear();
- _octree->GetRoot()->GetLeafsWith(Leafs,data->boundingbox[i],data->boundingbox[2+j],data->boundingbox[4+k]);
- if (Leafs.size()!=8){
- int min = Leafs[0]->BoxSize();
- for (unsigned int l=0;l<Leafs.size();l++){
- if (min>Leafs[l]->BoxSize()) min = Leafs[l]->BoxSize();
- }
- for (unsigned int l=0;l<Leafs.size();l++){
- if (min<=(Leafs[l]->BoxSize())/4){
- Leafs[l]->Divide();
- smoothed = false;
- }
- }
- }
- }
- }
- }
- return smoothed;
- }
-
- bool smoothed = true;
- for (int n = 0 ; n < 8 ; n++){
- if (!this->_children[n]->Smooth()) smoothed = false;
- }
- return smoothed;
-}
-
-// Octree research, box specialisation from root
-
-void Box::GetLeafsWith(std::vector<Box*> &Leafs, int x, int y, int z){
-
- if (!this->HaveChildren()){
- Box* parent = this;
- BoxData* data;
- while (parent->GetParent()!=NULL){
- data = parent->GetParent()->GetData();
- parent=parent->GetParent();
- }
- Leafs.push_back(this);
- return;
- }
-
-
- if ((x<=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y <= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z <= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
- _children[0]->GetLeafsWith(Leafs,x,y,z);
- }
-
- if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y <= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z <= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
- _children[1]->GetLeafsWith(Leafs,x,y,z);
- }
-
- if ((x<=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z <= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
- _children[2]->GetLeafsWith(Leafs,x,y,z);
- }
-
- if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z <= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
- _children[3]->GetLeafsWith(Leafs,x,y,z);
- }
-
- if ((x<=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y <= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z >= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
- _children[4]->GetLeafsWith(Leafs,x,y,z);
- }
-
- if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y <= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z >= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
- _children[5]->GetLeafsWith(Leafs,x,y,z);
- }
-
- if ((x<=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z >= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
- _children[6]->GetLeafsWith(Leafs,x,y,z);
- }
-
- if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) & (z >= (this->_data->boundingbox[5]+this->_data->boundingbox[4])/2) ){
- _children[7]->GetLeafsWith(Leafs,x,y,z);
- }
-
-}
diff --git a/contrib/arc/Classes/OctreeLSImage.h b/contrib/arc/Classes/OctreeLSImage.h
deleted file mode 100644
index f0835942841b938525f410e1bc4d9812f61d4608..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/OctreeLSImage.h
+++ /dev/null
@@ -1,128 +0,0 @@
-//
-// Description:
-//
-//
-// Author: <Boris Sedji>, 12/2009
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-
-#ifndef _OCTREELSIMAGE_H_
-#define _OCTREELSIMAGE_H_
-
-#include "itkImage.h"
-#include "itkImageFileReader.h"
-#include "itkImageFileWriter.h"
-#include "PView.h"
-#include "PViewData.h"
-
-
-//data defining a box
-
-class BoxData{
- public:
- int boundingbox[6]; // x1, x2, y1, y2, z1, z2
- float LevelSetValue[2][2][2]; // LevelSetValue[x][y][z]
-};
-
-class OctreeLSImage;
-class Box;
-
-// Octree class to decompose a float valued 3d image
-
-class OctreeLSImage
-{
-
- private :
-
- // Octree's root
- Box* _root;
- // Octree's image
- itk::Image< float, 3 >::Pointer _image;
- // Region size
- int _ImageSize[3];
- // Mesh list of nodes
- std::vector< std::vector <int> > _ListNodes;
- // Mesh list of elements
- std::vector< std::vector <int> > _ListElements;
- // Mesh list of levelsetvalue
- std::vector< float > _ListLSValue;
- // List of HangingNodes
- std::map< int ,std::vector<int> > _HangingNodes;
- // Number of box leafs
- int _LeafNumber;
- // Count octree's leafs
- void SetLeafNumber();
- void FillMeshInfo(int & pas);
- void FillHangingNodes();
-
- public :
-
- OctreeLSImage(itk::Image< float, 3 >::Pointer image);
- //~OctreeLSImage();
- void Erase();
- Box* GetRoot(){return _root;}
- // Refine with related mesh conditions (see implementation)
- void Mesh(int maxsize, int minsize);
- itk::Image< float, 3 >::Pointer GetImage(){return _image;}
- int GetLeafNumber(){this->SetLeafNumber();return _LeafNumber;}
- // Refine if too much generations between adjacent leafs
- bool Smooth();
- //std::vector< std::vector <int> > *GetListElements();
- // Create GModel
- GModel* CreateGModel(bool simplex, double facx, double facy, double facz,int & sizemax,int & sizemin);
- // Create PView representation of the level set
- PView* CreateLSPView(GModel* m);
- int* GetSize(){return _ImageSize;}
-
- std::vector< float >* getListLSValue(){return &_ListLSValue;}
- std::map< int ,std::vector<int> >* getListHangingNodes(){return &_HangingNodes;}
-};
-
-
-class Box{
- private :
-
- // children pointers
- Box* _children[8];
- Box* _parent;
- BoxData* _data;
- OctreeLSImage* _octree;
- int _ElementNodes[8]; // utilisé seulement par balayage
-
- public :
-
- Box(BoxData* data,Box* parent,OctreeLSImage* octree);
- BoxData* GetData(){return _data;}
- // Does it have children
- bool HaveChildren();
- // Fill level set value in data with image values
- void FillLevelSetValue(itk::Image< float, 3 >::Pointer image, BoxData *data);
- // The smallest length of the box
- int BoxSize();
- // Is the box crossed by the iso zero levelset
- bool IsHomogeneous();
- // Recursive function to refine until condition is reached (see implementation)
- void Mesh(int & maxsize, int & minsize);
- // Divide by eight the box leaf
- void Divide();
- void PrintLevelSetValue();
- Box* GetParent(){return _parent;};
- Box* GetChild(int n){return _children[n];};
- // Recursive function to count the leafs from the root
- void CountLeafNumber(int &LeafNumber);
- // Recursive function to smoothed the octree (see OctreeLSImage::Smooth)
- bool Smooth();
- // Give all the leafs containing this node
- void GetLeafsWith(std::vector<Box*> &Leafs, int x, int y, int z);
- // Recursive function to create the list of the mesh nodes, elements and levelset values
- void FillMeshInfo(std::vector<std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector< float > &ListLSValue,int &pas);
- // Recursive function to create the list of mesh elements in relation with the list of nodes
- void SetElementNode(int pos,int iti){_ElementNodes[pos]=iti;}
- void FillElementsNode(std::vector< std::vector<int> > &ListElements);
- void FillHangingNodes(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements,std::map< int ,std::vector<int> > &HangingNodes);
- };
-
-#endif
diff --git a/contrib/arc/Classes/OctreeSolver.cpp b/contrib/arc/Classes/OctreeSolver.cpp
deleted file mode 100644
index f34998b3b2ac40092077528e536836d739bc1b91..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/OctreeSolver.cpp
+++ /dev/null
@@ -1,486 +0,0 @@
-//
-// 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
deleted file mode 100644
index 4298f1e7c931221baac0442f179855e8f8ff913c..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/OctreeSolver.h
+++ /dev/null
@@ -1,53 +0,0 @@
-//
-// 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/QuadtreeLSImage.cpp b/contrib/arc/Classes/QuadtreeLSImage.cpp
deleted file mode 100644
index 4a49995bb2b81536905bad6ff3ce826f68b220ca..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/QuadtreeLSImage.cpp
+++ /dev/null
@@ -1,1017 +0,0 @@
-//
-// Description:
-//
-//
-// Author: <Boris Sedji>, 12/2009
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-#include "GModel.h"
-#include "QuadtreeLSImage.h"
-#include "MElement.h"
-
-QuadtreeLSImage::QuadtreeLSImage(itk::Image< float, 2 >::Pointer image){
-
- itk::Image< float, 2 >::RegionType region;
- // Define a region to get size of the image
- region = image->GetLargestPossibleRegion ();
-
- // Size will be a power of two larger than the image size
-
- int size[3];
- size[0] = 1;
- size[1] = 1;
-
- _ImageSize[0] = region.GetSize(0);
- _ImageSize[1] = region.GetSize(1);
-
- while (size[0]<region.GetSize(0)){
- size[0] = size[0]*2;
- }
- while (size[1]<region.GetSize(1)){
- size[1] = size[1]*2;
- }
-
- Box2DData *data = new Box2DData;
-
- // The window of the image root
-
- if (size[0]>= size[1])
- {
- // The window of the image root
- data->boundingbox[0]=0;
- data->boundingbox[1]=size[0];
- data->boundingbox[2]=0;
- data->boundingbox[3]=size[0];
-
- }
- if (size[1]>= size[0])
- {
- // The window of the image root
- data->boundingbox[0]=0;
- data->boundingbox[1]=size[1];
- data->boundingbox[2]=0;
- data->boundingbox[3]=size[1];
-
- }
-
- // initialisation of the quadtree
- _root = new Box2D(data,NULL,this);
- _root->FillLevelSetValue(image,data);
- _image = image;
- _LeafNumber = 0;
-
-}
-
-
-void QuadtreeLSImage::Mesh(int maxsize,int minsize){
-
- // taillemax //nombre de pixels max dans une direction
-
-
- if (this->_root->HaveChildren()){
- return;
- }
- if (_root->BoxSize()>maxsize || !_root->IsHomogeneous()){
- _root->Mesh( maxsize, minsize);
- };
-
-}
-
-
-Box2D::Box2D(Box2DData* data, Box2D* parent,QuadtreeLSImage* quadtree){
-
- _data = data;
- _parent = parent;
- _children[0] = NULL;
- _quadtree = quadtree;
-
-}
-
-bool Box2D::HaveChildren(){
-
- if (this->_children[0]==NULL) return false;
- else return true;
-
-}
-
-
-void Box2D::FillLevelSetValue(itk::Image< float, 2 >::Pointer image, Box2DData *data){
-
- float pixelValue;
- itk::Image< float, 2 >::IndexType pixelIndex;
- itk::Image< float, 2 >::RegionType region;
- region = image->GetLargestPossibleRegion ();
-
- // fill with pixel value, and if out of region size, fill with value at largest region size
- for (int i=0;i<2;i++){
- for (int j=0;j<2;j++){
-
- if (data->boundingbox[i] < region.GetSize(0)) pixelIndex[0] = data->boundingbox[i];
- else pixelIndex[0] = region.GetSize(0)-1; // x position
-
- if (data->boundingbox[2+j] < region.GetSize(1)) pixelIndex[1] = data->boundingbox[2+j];
- else pixelIndex[1] = region.GetSize(1)-1; // y position
-
- pixelValue = image->GetPixel( pixelIndex );
- data->LevelSetValue[i][j] = pixelValue;
-
- }
- }
-
-}
-
-int Box2D::BoxSize(){
-
- int SizeX;
- int SizeY;
- int size;
-
- SizeX = this->_data->boundingbox[1] - this->_data->boundingbox[0];
- SizeY = this->_data->boundingbox[3] - this->_data->boundingbox[2];
-
- if (SizeX<SizeY) size = SizeX;
- else size = SizeY;
-
- return size;
-
-}
-
-bool Box2D::IsHomogeneous(){
-
- float OneValue = this->_data->LevelSetValue[0][0];
-
- for (int i=0;i<2;i++){
- for (int j=0;j<2;j++){
- OneValue = OneValue*this->_data->LevelSetValue[i][j];
- if (OneValue<0) {
- return false; // if sign change => not homogeneous
- };
- OneValue = this->_data->LevelSetValue[i][j];
- }
- }
- return true;
-
-}
-
-
-// Recursive dividing function - eight nodes
-
-void Box2D::Mesh(int & maxsize, int & minsize){
-
-
- if (((this->BoxSize()<=maxsize) & (this->IsHomogeneous())) | (this->BoxSize()<=minsize)){ // (max size & homogenous) ou (min size) => dont divide
- return; // dont divide
- }
-
- // else we divide...
- this->Divide();
-
- for (int n = 0;n<4;n++){
- _children[n]->Mesh(maxsize,minsize);
- }
-
-}
-
-void Box2D::PrintLevelSetValue(){
-
- for (int i=0;i<2;i++){
- for (int j=0;j<2;j++){
- std::cout<<"\n "<< this->_data->LevelSetValue[i][j];
- }
- }
-
-}
-
-
-GModel *QuadtreeLSImage::CreateGModel(bool simplex, double facx, double facy,double facz, int & sizemax,int & sizemin){
-
- double eps = 1e-6;
- this->FillMeshInfo(sizemin);
- this->FillHangingNodes();
-
- std::cout<<"HangingNodes size :"<< _HangingNodes.size()<<"\n";
-
- std::map<int, MVertex*> vertexMap;
- std::vector<MVertex*> vertexVector;
-
- std::vector<int> numElement;
- std::vector<std::vector<int> > vertexIndices;
- std::vector<int> elementType ;
- std::vector<int> physical;
- std::vector<int> elementary;
- std::vector<int> partition;
- std::vector<double> d;
- std::map<int, std::vector<double> > data;
- data.clear();
-
- int numVertices;
-
- std::vector<std::vector<int> >::iterator ite;
-
- int i = 1;
- int k = 1;
-
- int a;
-
- int xlimit = _ImageSize[0] - _ImageSize[0]%sizemax; // lost of image pixels...
- int ylimit = _ImageSize[1] - _ImageSize[1]%sizemax;
-
- std::cout<<"\nxlimit :"<< xlimit<<"\n";
- std::cout<<"ylimit :"<< ylimit<<"\n";
-
- std::cout<<"Listnodes size :"<<_ListNodes.size()<<"\n";
-
- i = 1;
- ite = _ListElements.begin();
-
- if (!simplex) // first try, to be improved
- {
- while (ite!=_ListElements.end()) {
- int num, type, phys = 0, ele = 0, part = 0, numVertices;
- num = i;
- type = 3; // MSH_QUA_4
- ele = 1;
- phys = 1;
- part = 0;
- numVertices = MElement::getInfoMSH(type);
- std::vector <int> indices;
- for (int j = 0; j < numVertices; j++){
- indices.push_back((*ite)[j]);
- }
- numElement.push_back(i);
- vertexIndices.push_back(indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- indices.clear();
- ite++;
- i++;
- }
- }
- else // if simplex
- {
- while (ite!=_ListElements.end()) {
- bool inImage;
- inImage = true;
- for (int j = 0;j<4;j++)
- {
-// std::cout<<"_[(*ite)[j]]"<< (*ite)[j]<<"\n";
-// std::cout<<"_ListNodes[(*ite)[j]][0]"<< _ListNodes[(*ite)[j]-1][0] <<"\n";
-
- if ( _ListNodes[(*ite)[j]-1][0] > xlimit | _ListNodes[(*ite)[j]-1][1] > ylimit )
- {
- inImage = false;
- }
- }
-
- //inImage = true;
- if (inImage)
- {
-
- k++;
- for (int j = 0;j<4;j++)
- {
- if (vertexMap.find((*ite)[j]) == vertexMap.end()) // if not in
- {
- float xyz[3];
- d.clear();
- xyz[0] = (_ListNodes[(*ite)[j]-1][0])*facx;
- xyz[1] = (_ListNodes[(*ite)[j]-1][1])*facy;
- xyz[2] = 0;
- MVertex* newVertex = new MVertex(xyz[0], xyz[1], xyz[2], 0, (*ite)[j]);
- vertexMap[(*ite)[j]] = newVertex;
- d.push_back(_ListLSValue[(*ite)[j]-1]);
- data[(*ite)[j]]=d;
- d.clear();
- }
- }
-
- int num, type, phys = 0, ele = 0, part = 0, numVertices;
- num = i;
- numVertices = 4;
- std::vector <int> indices;
- std::vector<int> front_indices;
-
- // --- domain elements - first triangle ---
-
-
- indices.clear();
- for (int j = 0; j < 3; j++){
- indices.push_back((*ite)[j]);
- }
- type = 2; // MSH_TRI_3
- ele = 5;
- phys = 5;
- part = 0;
- numElement.push_back(i);
- vertexIndices.push_back(indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
-
- // --- frontier element face 1 -> y = 0 ---
-
- front_indices.clear();
- for (int j =0;j<3;j++)
- {
- if (vertexMap[indices[j]]->y() < eps*facy)
- {
- front_indices.push_back(indices[j]);
- }
- }
- if (front_indices.size() == 2)
- {
- type = 1; // MSH_LINE
- ele = 1;
- phys = 1;
- part = 0;
- numVertices = 2;
- numElement.push_back(i);
- vertexIndices.push_back(front_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
- // --- frontier element face 2 -> x = 0 ---
-
- front_indices.clear();
- for (int j =0;j<3;j++)
- {
- if (vertexMap[indices[j]]->x() < eps*facx)
- {
- front_indices.push_back(indices[j]);
- }
- }
- if (front_indices.size() == 2)
- {
- type = 1; // MSH_LINE
- ele = 2;
- phys = 2;
- part = 0;
- numVertices = 2;
- numElement.push_back(i);
- vertexIndices.push_back(front_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
-
-
- // - frontier element face 3 -> y = _size[1]
- front_indices.clear();
- for (int j =0;j<3;j++)
- {
- if (vertexMap[indices[j]]->y() > (ylimit - eps)*facy)
- {
- front_indices.push_back(indices[j]);
- }
- }
- if (front_indices.size() == 2)
- {
- type = 1; // MSH_LINE
- ele = 3;
- phys = 3;
- part = 0;
- numVertices = 2;
- numElement.push_back(i);
- vertexIndices.push_back(front_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
- // --- frontier element face 4 -> x = _size[0] ---
-
- front_indices.clear();
- for (int j =0;j<3;j++)
- {
- if (vertexMap[indices[j]]->x() > (xlimit-eps)*facx)
- {
- front_indices.push_back(indices[j]);
- }
- }
- if (front_indices.size() == 2)
- {
- type = 1; // MSH_LINE
- ele = 4;
- phys = 4;
- part = 0;
- numVertices = 2;
- numElement.push_back(i);
- vertexIndices.push_back(front_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
-
-
- indices.clear();
- // second triangle
- for (int j = 2; j < 5; j++){
- indices.push_back((*ite)[j%4]);
- }
- type = 2; // MSH_TRI_3
- ele = 5;
- phys = 5;
- part = 0;
- numElement.push_back(i);
- vertexIndices.push_back(indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
-
- // --- frontier element face 1 -> y = 0 ---
- front_indices.clear();
- for (int j =0;j<3;j++)
- {
- if (vertexMap[indices[j]]->y() < eps*facy)
- {
- front_indices.push_back(indices[j]);
- }
- }
- if (front_indices.size() == 2)
- {
- type = 1; // MSH_LINE
- ele = 1;
- phys = 1;
- part = 0;
- numVertices = 2;
- numElement.push_back(i);
- vertexIndices.push_back(front_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
-
- // --- frontier element face 2 -> x = 0 ---
-
- front_indices.clear();
- for (int j =0;j<3;j++)
- {
- if (vertexMap[indices[j]]->x() < eps*facx)
- {
- front_indices.push_back(indices[j]);
- }
- }
- if (front_indices.size() == 2)
- {
- type = 1; // MSH_LINE
- ele = 2;
- phys = 2;
- part = 0;
- numVertices = 2;
- numElement.push_back(i);
- vertexIndices.push_back(front_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
-
- // --- frontier element face 3 -> y = _size[1] ---
- front_indices.clear();
- for (int j =0;j<3;j++)
- {
- if (vertexMap[indices[j]]->y() > (ylimit - eps)*facy)
- {
- front_indices.push_back(indices[j]);
- }
- }
- if (front_indices.size() == 2)
- {
- type = 1; // MSH_LINE
- ele = 3;
- phys = 3;
- part = 0;
- numVertices = 2;
- numElement.push_back(i);
- vertexIndices.push_back(front_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
- front_indices.clear();
-
-
- // --- frontier element face 4 -> x = _size[0] ---
-
- front_indices.clear();
- for (int j =0;j<3;j++)
- {
- if (vertexMap[indices[j]]->x() > (xlimit-eps)*facx)
- {
- front_indices.push_back(indices[j]);
- }
- }
- if (front_indices.size() == 2)
- {
- type = 1; // MSH_LINE
- ele = 4;
- phys = 4;
- part = 0;
- numVertices = 2;
- numElement.push_back(i);
- vertexIndices.push_back(front_indices);
- elementType.push_back(type);
- physical.push_back(phys);
- elementary.push_back(ele);
- partition.push_back(part);
- i++;
- }
-
- }
- ite++; // next quadtree leaf
- }
- }
-
-
- std::cout<<"numElement size :"<<numElement.size()<<"\n";
- std::cout<<"vertexIndices size :"<<vertexIndices.size()<<"\n";
-
-
- std::cout<<"data size :"<<data.size()<<"\n";
- std::cout<<"vertexMap size :"<<vertexMap.size()<<"\n";
- GModel *gmod = GModel::createGModel(vertexMap,numElement,vertexIndices,elementType,physical,elementary,partition);
- // Write a .msh file with the level set values as postview data
- std::string postTypeName = "LevelSet Value" ;
- PView *pv = new PView (postTypeName, "NodeData", gmod, data);
- //PView *pv = octree.CreateLSPView(m);
- bool useadapt = true;
- pv->getData(useadapt)->writeMSH("LSPView.msh", false);
-
- std::cout<<"getNumScalars :"<< pv->getData(useadapt)->getNumScalars()<<"\n";
-
- return gmod;
-
-}
-
-PView *QuadtreeLSImage::CreateLSPView(GModel* m){
-
- std::map<int, std::vector<double> > data;
- std::vector<float>::iterator itf;
- itf = _ListLSValue.begin();
- int i = 1;
-
- while (itf!=_ListLSValue.end()){
- std::vector<double> d;
- d.push_back((*itf)) ;
- data[i] = d;
- d.clear();
- itf++;
- i++;
- }
- std::string postTypeName = "LevelSet Value" ;
- PView *pv = new PView (postTypeName, "NodeData", m, data, 0.0);
-
- return pv;
-
-}
-
-
-
-void QuadtreeLSImage::FillHangingNodes(){
-
- _root->FillHangingNodes(_ListNodes,_ListElements,_HangingNodes);
-
-}
-
-
-
-
-
-void Box2D::FillHangingNodes(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements,std::map< int ,std::vector<int> > &HangingNodes){
-
-
- std::vector<std::vector<int> >::iterator it;
-
- it = ListNodes.begin();
- int k=0;
- while (it!=ListNodes.end())
- {
- std::vector<Box2D*> Leafs;
- std::vector<int> MasterNodes;
- Leafs.clear();
- MasterNodes.clear();
- GetLeafsWith(Leafs, (*it)[0], (*it)[1]);
- if (Leafs.size()==3)
- {
- Box2D *MaxLeaf;
- MaxLeaf = Leafs[0];
- for (int i=1;i<3;i++)
- {
- if (MaxLeaf->BoxSize() < Leafs[i]->BoxSize())
- MaxLeaf = Leafs[i];
- }
- for (int i = 0 ; i<4;i++)
- {
- if ( ListNodes[MaxLeaf->_ElementNodes[i]-1][0] == (*it)[0] | ListNodes[MaxLeaf->_ElementNodes[i]-1][1]== (*it)[1])
- {
- MasterNodes.push_back(MaxLeaf->_ElementNodes[i]);
- }
- }
- HangingNodes[k+1] = MasterNodes;
- }
- it++;
- k++;
- }
-
-}
-
-
-
-
-void QuadtreeLSImage::FillMeshInfo(int & pas){
-
- if (!_ListNodes.empty()){
- return;
- }
-
- _root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue,pas);
-
-}
-
-
-
-// fonction par balayage de région
-
-
-void Box2D::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector<float> &ListLSValue,int &pas){
-
-
- int size = _quadtree->GetRoot()->BoxSize();
-
- int xpas =pas; // Ã modifier si sizeZ n'est pas le plus petit
- int ypas = pas;
-
- std::vector<std::vector<int> >::iterator it;
-
- it = ListNodes.begin();
- int iti=0;
-
-
- for (int y = 0;y<=size;y+=ypas)
- for (int x = 0;x<=size; x+=xpas)
- {
- std::vector<Box2D*> Leafs;
- Leafs.clear();
- GetLeafsWith(Leafs, x, y);
- std::map<int,std::vector< int > > ElementsNodes;
- std::vector<int> NodesIn;
-
- bool added;
- added = false;
- bool HangingNodes;
- HangingNodes = false;
- for (unsigned int l=0;l<Leafs.size();l++){
- if (!added){
- std::vector<int> XYZ;
- XYZ.push_back(x);
- XYZ.push_back(y);
- XYZ.push_back(0);
- for (int i = 0 ; i<2;i++ )
- for (int j = 0 ; j<2;j++ )
- if (x == Leafs[l]->GetData()->boundingbox[i] & y == Leafs[l]->GetData()->boundingbox[j+2] ){
- ListLSValue.push_back(Leafs[l]->GetData()->LevelSetValue[i][j]);
- ListNodes.push_back(XYZ);
- added = true;
- iti++;
- }
-
- }else break;
- }
- if (added){
- for (unsigned int l=0;l<Leafs.size();l++){
- for (int i = 0 ; i<2;i++ )
- for (int j = 0 ; j<2;j++ ){
- int pos ;
- pos = i*2+j*1;
- if (x == Leafs[l]->GetData()->boundingbox[i] & y == Leafs[l]->GetData()->boundingbox[j+2]){
- Leafs[l]->SetElementNode(pos,iti);
- }
- }
- }
- }
- }
- std::cout<<"\nNumber Nodes : "<< ListNodes.size();
- _quadtree->GetRoot()->FillElementsNode(ListElements);
-
-}
-
-void Box2D::FillElementsNode(std::vector< std::vector<int> > &ListElements){
-
- // if it s a leaf then fill list node if node dont exist in list
- if (!this->HaveChildren()) {
- std::vector<int> ElementNodes;
- int temp;
- temp = _ElementNodes[3];
- _ElementNodes[3] = _ElementNodes[2];
- _ElementNodes[2] = temp;
- for (int i = 0 ; i < 4 ; i++) ElementNodes.push_back(_ElementNodes[i]);
- ListElements.push_back(ElementNodes);
-// std::cout<<"\n Nodes :";
-// for (int i = 0 ; i < 4 ; i++)
-// std::cout<<" " << ElementNodes[i];
- return;
- }
-
- for (int n = 0 ; n < 4 ; n++){
- this->_children[n]->FillElementsNode(ListElements);
- }
-
-}
-
-
-
-
-/*
-void Box::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector<float> &ListLSValue){
-
- // if it s a leaf then fill list node if node dont exist in list
- if (!this->HaveChildren()) {
- std::vector<int> NodesIn;
- NodesIn.clear();
- int temp;
- for (int i = 0 ; i<2;i++ ){
- for (int k = 0 ; k<2;k++ )
- for (int j = 0 ; j<2;j++ ){
- std::vector<int> XYZ;
- float LSValue;
- XYZ.push_back(this->GetData()->boundingbox[i]);
- XYZ.push_back(this->GetData()->boundingbox[2+j]);
- XYZ.push_back(this->GetData()->boundingbox[4+k]);
- LSValue = this->GetData()->LevelSetValue[i][j][k];
- std::vector<std::vector<int> >::iterator it;
- it = ListNodes.begin();
- bool add = true;
- int l;
- l = 1;
- while (it!=ListNodes.end()){ // attention segmentation & pas opti
- if (((*it)[0] == XYZ[0]) & ((*it)[1] == XYZ[1]) & ((*it)[2] == XYZ[2])){
- add = false;
- NodesIn.push_back(l);
- break;
- }
- it++;
- l++;
- }
- if (add) {
- ListNodes.push_back(XYZ);
- ListLSValue.push_back(LSValue);
- NodesIn.push_back(l);
- }
- XYZ.clear();
- }
- }
- temp = NodesIn[3];
- NodesIn[3] = NodesIn[2];
- NodesIn[2] = temp;
- temp = NodesIn[7];
- NodesIn[7] = NodesIn[6];
- NodesIn[6] = temp;
- ListElements.push_back(NodesIn);
- return;
- }
-
- for (int n = 0 ; n < 8 ; n++){
- this->_children[n]->FillMeshInfo(ListNodes,ListElements,ListLSValue);
- }
-
-}*/
-
-/*
-std::vector< std::vector <int> > *QuadtreeLSImage::GetListElements(){
-
- if (!_ListElements.empty()){
- return &_ListElements;
- }
-
- _root->GetLeafElements(_ListElements);
-
- return &_ListElements;
-
-}*/
-
-
-/*
-void Box2D::GetLeafElements(std::vector< std::vector<int> > &ListElements){
-
-// if (!this->HaveChildren()) {
-// std::vector<std::vector<int> >::iterator it;
-// std::vector<std::vector<int> > *ListNodes;
-// ListNodes = this->_quadtree->GetListNodes();
-// it = ListNodes->begin();
-// std::vector<int> NodesIn;
-// int l = 1;
-// while (it!=ListNodes->end()) {
-// for (int i = 0 ; i<2;i++ ){
-// for (int j = 0 ; j<2;j++ ){
-// for (int k = 0 ; k<2;k++ ){
-// if ((_data->boundingbox[i] == (*it)[0]) & (_data->boundingbox[2+j] == (*it)[1]) & (_data->boundingbox[4+k] == (*it)[2]))
-// NodesIn.push_back(l);
-// }
-// }
-// }
-// it++;
-// l++;
-// }
-// ListElements.push_back(NodesIn);
-// return;
-// }
-//
-// for (int n = 0 ; n < 8 ; n++){
-// this->_children[n]->GetLeafElements(ListElements);
-// }
-
-}
-*/
-
-
-
-
-void QuadtreeLSImage::SetLeafNumber(){
- _LeafNumber = 0;
- _root->CountLeafNumber(_LeafNumber);
-}
-
-void Box2D::CountLeafNumber(int &LeafNumber){
-
- if (this->GetChild(0)==NULL) {
- LeafNumber = LeafNumber + 1;
- return;
- }
-
- for (int n = 0 ; n < 4 ; n++){
- this->_children[n]->CountLeafNumber(LeafNumber);
- }
-
-}
-
-// Divide function
-
-void Box2D::Divide(){
-
-// box divide in 4 boxes with limits xl = (xmin + xmax)/2 ; yl = (ymin + ymax)/2
-
- int n;
-
- // Box 1 : < xl ; < yl
-
- n = 0;
- Box2DData* data1 = new Box2DData;
-
- data1->boundingbox[0]=this->_data->boundingbox[0];
- data1->boundingbox[1]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data1->boundingbox[2]=this->_data->boundingbox[2];
- data1->boundingbox[3]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
-
-
- this->_children[n] = new Box2D(data1,this,this->_quadtree);
- this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
-
- // Box 2 : > xl ; < yl
-
- n = 1;
- Box2DData* data2 = new Box2DData;
-
- data2->boundingbox[0]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data2->boundingbox[1]=this->_data->boundingbox[1];
- data2->boundingbox[2]=this->_data->boundingbox[2];
- data2->boundingbox[3]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
-
-
- this->_children[n] = new Box2D(data2,this,this->_quadtree);
- this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
-
- // Box 3 : < xl ; > yl
-
- n = 2;
- Box2DData* data3 = new Box2DData;
-
- data3->boundingbox[0]=this->_data->boundingbox[0];
- data3->boundingbox[1]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data3->boundingbox[2]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
- data3->boundingbox[3]=this->_data->boundingbox[3];
-
-
- this->_children[n] = new Box2D(data3,this,this->_quadtree);
- this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
-
- // Box 4 : > xl ; > yl
-
- n = 3;
- Box2DData* data4 = new Box2DData;
-
- data4->boundingbox[0]=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2;
- data4->boundingbox[1]=this->_data->boundingbox[1];
- data4->boundingbox[2]=(this->_data->boundingbox[3]+this->_data->boundingbox[2])/2;
- data4->boundingbox[3]=this->_data->boundingbox[3];
-
-
- this->_children[n] = new Box2D(data4,this,this->_quadtree);
- this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
-
-
-}
-
-bool QuadtreeLSImage::Smooth(){
-
- return _root->Smooth();
-
-}
-
-
-bool Box2D::Smooth(){
-
-
- if (!this->HaveChildren()) {
- bool smoothed = true;
- Box2DData*data = this->GetData();
- std::vector<Box2D*> Leafs;
- for (int i = 0 ; i<2;i++ ){
- for (int j = 0 ; j<2;j++ ){
- Leafs.clear();
- _quadtree->GetRoot()->GetLeafsWith(Leafs,data->boundingbox[i],data->boundingbox[2+j]);
- if (Leafs.size()!=4){
- int min = Leafs[0]->BoxSize();
- for (unsigned int l=0;l<Leafs.size();l++){
- if (min>Leafs[l]->BoxSize()) min = Leafs[l]->BoxSize();
- }
- for (unsigned int l=0;l<Leafs.size();l++){
- if (min<=(Leafs[l]->BoxSize())/4){
- Leafs[l]->Divide();
- smoothed = false;
- }
- }
- }
- }
- }
- return smoothed;
- }
-
- bool smoothed = true;
- for (int n = 0 ; n < 4 ; n++){
- if (!this->_children[n]->Smooth()) smoothed = false;
- }
- return smoothed;
-}
-
-// Quadtree research, box specialisation from root
-
-void Box2D::GetLeafsWith(std::vector<Box2D*> &Leafs, int x, int y){
-
- if (!this->HaveChildren()){
- Box2D* parent = this;
- Box2DData* data;
- while (parent->GetParent()!=NULL){
- data = parent->GetParent()->GetData();
- parent=parent->GetParent();
- }
- Leafs.push_back(this);
- return;
- }
-
-
- if ((x<=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y <= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) ){
- _children[0]->GetLeafsWith(Leafs,x,y);
- }
-
- if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y <= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) ){
- _children[1]->GetLeafsWith(Leafs,x,y);
- }
-
- if ((x<=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2) ){
- _children[2]->GetLeafsWith(Leafs,x,y);
- }
-
- if ((x>=(this->_data->boundingbox[1]+this->_data->boundingbox[0])/2) & (y >= (this->_data->boundingbox[3]+this->_data->boundingbox[2])/2)){
- _children[3]->GetLeafsWith(Leafs,x,y);
- }
-
-
-}
diff --git a/contrib/arc/Classes/QuadtreeLSImage.h b/contrib/arc/Classes/QuadtreeLSImage.h
deleted file mode 100644
index d7f3eb0257ac3a8d42de4a49cc13315ed65ae933..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/QuadtreeLSImage.h
+++ /dev/null
@@ -1,129 +0,0 @@
-//
-// Description:
-//
-//
-// Author: <Boris Sedji>, 12/2009
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-
-#ifndef _QUADTREELSIMAGE_H_
-#define _QUADTREELSIMAGE_H_
-
-#include "itkImage.h"
-#include "itkImageFileReader.h"
-#include "itkImageFileWriter.h"
-#include "PView.h"
-#include "PViewData.h"
-
-
-//data defining a box
-
-class Box2DData{
- public:
- int boundingbox[4]; // x1, x2, y1, y2
- float LevelSetValue[2][2]; // LevelSetValue[x][y][z]
-};
-
-class QuadtreeLSImage;
-class Box2D;
-
-// Quadtree class to decompose a float valued 3d image
-
-class QuadtreeLSImage
-{
-
- private :
-
- // Quadtree's root
- Box2D* _root;
- // Quadtree's image
- itk::Image< float, 2 >::Pointer _image;
- // Region size
- int _ImageSize[2];
- // Mesh list of nodes
- std::vector< std::vector <int> > _ListNodes;
- // Mesh list of elements
- std::vector< std::vector <int> > _ListElements;
- // Mesh list of levelsetvalue
- std::vector< float > _ListLSValue;
- // List of HangingNodes
- std::map< int ,std::vector<int> > _HangingNodes;
- // Number of box leafs
- int _LeafNumber;
- // Count octree's leafs
- void SetLeafNumber();
- void FillMeshInfo(int & pas);
- void FillHangingNodes();
-
- public :
-
- QuadtreeLSImage(itk::Image< float, 2 >::Pointer image);
- //~QuadtreeLSImage();
- void Erase();
- Box2D* GetRoot(){return _root;}
- // Refine with related mesh conditions (see implementation)
- void Mesh(int maxsize, int minsize);
- itk::Image< float, 2 >::Pointer GetImage(){return _image;}
- int GetLeafNumber(){this->SetLeafNumber();return _LeafNumber;}
- // Refine if too much generations between adjacent leafs
- bool Smooth();
- //std::vector< std::vector <int> > *GetListElements();
- // Create GModel
- GModel* CreateGModel(bool simplex, double facx, double facy,double facz,int & sizemax,int & sizemin);
- // Create PView representation of the level set
- PView* CreateLSPView(GModel* m);
- int* GetSize(){return _ImageSize;}
-
- std::vector< float >* getListLSValue(){return &_ListLSValue;}
- std::map< int ,std::vector<int> >* getListHangingNodes(){return &_HangingNodes;}
-
-};
-
-
-class Box2D{
- private :
-
- // children pointers
- Box2D* _children[4];
- Box2D* _parent;
- Box2DData* _data;
- QuadtreeLSImage* _quadtree;
- int _ElementNodes[4]; // utilisé seulement par balayage
-
- public :
-
- Box2D(Box2DData* data,Box2D* parent,QuadtreeLSImage* octree);
- Box2DData* GetData(){return _data;}
- // Does it have children
- bool HaveChildren();
- // Fill level set value in data with image values
- void FillLevelSetValue(itk::Image< float, 2 >::Pointer image, Box2DData *data);
- // The smallest length of the box
- int BoxSize();
- // Is the box crossed by the iso zero levelset
- bool IsHomogeneous();
- // Recursive function to refine until condition is reached (see implementation)
- void Mesh(int & maxsize, int & minsize);
- // Divide by eight the box leaf
- void Divide();
- void PrintLevelSetValue();
- Box2D* GetParent(){return _parent;};
- Box2D* GetChild(int n){return _children[n];};
- // Recursive function to count the leafs from the root
- void CountLeafNumber(int &LeafNumber);
- // Recursive function to smoothed the octree (see QuadtreeLSImage::Smooth)
- bool Smooth();
- // Give all the leafs containing this node
- void GetLeafsWith(std::vector<Box2D*> &Leafs, int x, int y);
- // Recursive function to create the list of the mesh nodes, elements and levelset values
- void FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector<float> &ListLSValue,int &pas);
- // Recursive function to create the list of mesh elements in relation with the list of nodes
- void SetElementNode(int pos,int iti){_ElementNodes[pos]=iti;}
- void FillElementsNode(std::vector< std::vector<int> > &ListElements);
- void FillHangingNodes(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements,std::map< int ,std::vector<int> > &HangingNodes);
- };
-
-#endif
diff --git a/contrib/arc/Classes/VectorXFEMFS.cpp b/contrib/arc/Classes/VectorXFEMFS.cpp
deleted file mode 100644
index 268f8f26d1a5f7ef67e86a2b644aa23f83de2884..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/VectorXFEMFS.cpp
+++ /dev/null
@@ -1,198 +0,0 @@
-//
-// 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)
-{
- std::set<int>::iterator it;
- it = _TagEnrichedVertex->find(ver->getNum());
- if (it==_TagEnrichedVertex->end()) 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++)
- {
- std::set<int>::iterator it;
- it = _TagEnrichedVertex->find(ele->getVertex(i)->getNum());
- if (it!=_TagEnrichedVertex->end()) 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++)
- {
- std::set<int>::iterator it;
- it = _TagEnrichedVertex->find(ele->getVertex(i)->getNum());
- if (it!=_TagEnrichedVertex->end())
- {
- 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++)
- {
- std::set<int>::iterator it;
- it = _TagEnrichedVertex->find(ele->getVertex(i)->getNum());
- if (it!=_TagEnrichedVertex->end())
- {
- 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++)
- {
- std::set<int>::iterator it;
- it = _TagEnrichedVertex->find(ele->getVertex(i)->getNum());
- if (it!=_TagEnrichedVertex->end()) // 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++)
- {
- std::set<int>::iterator it;
- it = _TagEnrichedVertex->find(ele->getVertex(i)->getNum());
- if (it!=_TagEnrichedVertex->end()) ndofs = ndofs + 1;
- }
- }
- return ndofs;
-}
-
-int ScalarLagrangeToXfemFS::getNumKeys(MVertex *ver)
-{
- std::set<int>::iterator it;
- it = _TagEnrichedVertex->find(ver->getNum());
- if (it!=_TagEnrichedVertex->end()) return 2 ;
- else return 1; // 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++)
- {
- std::set<int>::iterator it;
- it = _TagEnrichedVertex->find(ele->getVertex(i)->getNum());
- if (it!=_TagEnrichedVertex->end()) 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)
-{
- std::set<int>::iterator it;
- it = _TagEnrichedVertex->find(ver->getNum());
- if (it!=_TagEnrichedVertex->end())
- {
- 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
deleted file mode 100644
index 8a156bb5c88c7ff64a7ca897249874da3b5bded4..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/VectorXFEMFS.h
+++ /dev/null
@@ -1,93 +0,0 @@
-//
-// 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 <set>
-
-
-#include "functionSpace.h"
-#include "simpleFunction.h"
-
-
-class ScalarLagrangeToXfemFS : public ScalarLagrangeFunctionSpace{
- public:
-
- typedef TensorialTraits<double>::ValType ValType;
- typedef TensorialTraits<double>::GradType GradType;
-
- protected:
-
- std::set<int > *_TagEnrichedVertex;
- simpleFunction<double> *_funcEnrichment;
-
- public:
-//
- ScalarLagrangeToXfemFS(int i, std::set<int > &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;
-
-
- ScalarXFEMToVectorFS(int id , std::set<int > & 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::set<int > &TagEnrichedVertex , simpleFunction<double> * funcEnrichment) :
- ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment),
- BasisVectors[comp1],comp1) {}
-
- ScalarXFEMToVectorFS(int id,Along comp1,Along comp2, std::set<int > &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::set<int > &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
deleted file mode 100644
index c4058129261ba8f5e13fbaa034b6344fd6d171ee..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/XFEMelasticitySolver.cpp
+++ /dev/null
@@ -1,191 +0,0 @@
-//
-// 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"
-#include "MPoint.h"
-#endif
-
-
-#include "DILevelset.h"
-//#include "VectorXFEMFS.cpp"
-#include "XFEMelasticitySolver.h"
-#include "FuncHeaviside.h"
-
-#include "NodeEnrichedFS.cpp"
-
-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.insert(e->getParent()->getVertex(k)->getNum());
- }
- }
- }
- }
-
- // 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
- 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 = 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
-
- 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);
- for (int i=0;i<pAssembler->sizeOfR();i++) std::cout<< lsys->getFromSolution(i) << "\n";
-
-}
-
-
-
-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;
- 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/XFEMelasticitySolver.h b/contrib/arc/Classes/XFEMelasticitySolver.h
deleted file mode 100644
index 57266017f95091e9e0466d55c1035844b2df5822..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/XFEMelasticitySolver.h
+++ /dev/null
@@ -1,38 +0,0 @@
-//
-// 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::set<int > _TagEnrichedVertex;
- // enriched comp
- std::vector<int> _EnrichComp;
- // 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/Classes/gLevelSetMesh.cpp b/contrib/arc/Classes/gLevelSetMesh.cpp
deleted file mode 100644
index c83128397140a328a0e3a8787347b8a953b22ace..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/gLevelSetMesh.cpp
+++ /dev/null
@@ -1,21 +0,0 @@
-//
-// 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
deleted file mode 100644
index 57ed49faf0790fdcd72bded4cc22bb6dec451e0a..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/gLevelSetMesh.h
+++ /dev/null
@@ -1,61 +0,0 @@
-//
-// 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
deleted file mode 100644
index 06fc20fadcc9e2a2350241e047771d64d6fe53f1..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/xFemFunctionSpace.cpp
+++ /dev/null
@@ -1,220 +0,0 @@
-//
-// 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
deleted file mode 100644
index b081870d5aed5452b7975e3cac99fa541f69bdc7..0000000000000000000000000000000000000000
--- a/contrib/arc/Classes/xFemFunctionSpace.h
+++ /dev/null
@@ -1,111 +0,0 @@
-//
-// 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
diff --git a/contrib/arc/ImportLS2dImage.cpp b/contrib/arc/ImportLS2dImage.cpp
deleted file mode 100755
index 1ac56bda2caf203b5a378d384422d31ab4ccc6f9..0000000000000000000000000000000000000000
--- a/contrib/arc/ImportLS2dImage.cpp
+++ /dev/null
@@ -1,87 +0,0 @@
-//
-// 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++;
- }
-
- bool simplex = 1;
- double facx = 1;
- double facy = 1;
- // Create GModel with the octree mesh
- GModel *m = quadtree.CreateGModel(simplex,facx,facy,sizemax, sizemin);
-
- // 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/ImportLSImage.cpp b/contrib/arc/ImportLSImage.cpp
deleted file mode 100644
index 715a80504daa465466edcb3f946a0d887580c09e..0000000000000000000000000000000000000000
--- a/contrib/arc/ImportLSImage.cpp
+++ /dev/null
@@ -1,92 +0,0 @@
-//
-// Description:
-//
-//
-// Author: <Boris Sedji>, 12/2009
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-
-#ifndef _ITK_H_
-
-#include "itkImage.h"
-#include "itkImageFileReader.h"
-#include "itkImageFileWriter.h"
-#define _ITK_H_
-
-#endif
-
-#include "Octree.h"
-#include <stdio.h>
-#include "Gmsh.h"
-#include "GModel.h"
-#include "MElement.h"
-#include "gmshLevelset.h"
-#include "OctreeLSImage/OctreeLSImage.cpp"
-
-
-
-int main( int argc, char *argv[] )
-{
- if( argc < 4 )
- {
- std::cerr << "Missing Parameters " << std::endl;
- std::cerr << "Usage: " << argv[0];
- std::cerr << " inputLevelSetImage MeshSizeMax MeshSizeMin"<< std::endl;
- return 1;
- }
-
- typedef float PixelTypeFloat;
- const unsigned int Dimension = 3;
- 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) << " x " << region.GetSize(2);
-
- OctreeLSImage octree(image);
-
- int sizemax = atoi(argv[2]);
- int sizemin = atoi(argv[3]);
-
- // Mesh with conditions on mesh size
- octree.Mesh(sizemax,sizemin);
-
- // Smoothing mesh to avoid too much generation between adjacent elements
- bool statut=false;
- int k = 0;
- std::cout<<"\nLeaf Number : "<< (octree.GetLeafNumber())<<"\n";
- while((!statut) & (k < 20)){
- statut = octree.Smooth();
- std::cout<<"\nk : "<< k;
- std::cout<<"\nsmoothed : " << (int)statut;
- std::cout<<"\nLeaf Number : "<< (octree.GetLeafNumber())<<"\n";
- k++;
- }
-
- // Create GModel with the octree mesh
- GModel *m = octree.CreateGModel();
-
- // Write a .msh file with the mesh
- std::string ModelName = "OctreeMesh.msh" ;
- m->writeMSH(ModelName,2.1,false,false);
-
- // Write a .msh file with the level set values as postview data
- PView *pv = octree.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
deleted file mode 100644
index e01cd7289e7a411631c26b6d0b8464d1931a2435..0000000000000000000000000000000000000000
--- a/contrib/arc/XFEMElasticity.cpp
+++ /dev/null
@@ -1,55 +0,0 @@
-//
-// 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\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();
-
-}
-
diff --git a/contrib/arc/XFEMInclusion.cpp b/contrib/arc/XFEMInclusion.cpp
deleted file mode 100644
index 63ce7b37827ca259377bda86c78cce1d4596e635..0000000000000000000000000000000000000000
--- a/contrib/arc/XFEMInclusion.cpp
+++ /dev/null
@@ -1,114 +0,0 @@
-//
-// Description :
-//
-//
-// Author: <Boris Sedji>, 01/2010
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-
-#include "GModel.h"
-#include "Gmsh.h"
-#include "elasticitySolver.h"
-#include "PView.h"
-#include "PViewData.h"
-#include "highlevel.h"
-#include "groupOfElements.h"
-#include <iterator>
-
-
-#include "Classes/OctreeSolver.cpp"
-
-
-
-int main (int argc, char* argv[])
-{
-
-
-//
- if (argc < 2)
- {
- printf("usage : elasticity input_file_name.msh \n");
- return -1;
- }
-
- // globals are still present in Gmsh
- GmshInitialize(argc, argv);
-
- OctreeSolver mySolver(1000);
-
- 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();
-
-//
-//
-// // what is the input file ?
-//
-// if (argc != 2){
-// printf("Input file missing...\n");
-// return -1;
-// }
-//
-// // globals are still present in Gmsh
-//
-// GmshInitialize(argc, argv);
-//
-// // declare a GModel
-//
-// GModel *m = new GModel();
-//
-// // read from input file
-//
-// m->readMSH(argv[1]);
-//
-// // number of elements
-//
-// int nbr_Elements;
-// nbr_Elements = m->getNumMeshElements();
-// std::cout<< "\nLe nombre d'éléments du modèle non coupé est " << nbr_Elements;
-//
-// // level set plane definition
-//
-// double a(0.), b(1.), c(0.), d(-300); // normal vector and distance
-// int n(1); // tag
-//
-// gLevelset *ls = new gLevelsetPlane(a, b, c, d, n); // reference argument
-//
-// // cut model
-//
-// GModel *mm1 = m->buildCutGModel(ls);
-//
-// // save cut model in a file by the name of uncut + "_cut"
-//
-// std::string ModelName = argv[1] ;
-// ModelName.erase(ModelName.find(".",0),4);
-// ModelName.insert(ModelName.length(),"_cut.msh");
-// mm1->writeMSH(ModelName,2.1,false,false);
-//
-//
-// delete m;
-//
-// // stop gmsh
-//
-// GmshFinalize();
-
-//
-
-
-
-
-}
-
diff --git a/contrib/arc/boris.cpp b/contrib/arc/boris.cpp
deleted file mode 100644
index 5c7f62b6013dcde3cdd402296ccf7ad0915a0577..0000000000000000000000000000000000000000
--- a/contrib/arc/boris.cpp
+++ /dev/null
@@ -1,66 +0,0 @@
-
-#include "itkCurvatureAnisotropicDiffusionImageFilter.h"
-
-#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"
-#include "itkSigmoidImageFilter.h"
-
-#include "itkImage.h"
-#include "itkFastMarchingImageFilter.h"
-
-#include "itkBinaryThresholdImageFilter.h"
-
-#include "itkImageFileReader.h"
-#include "itkImageFileWriter.h"
-
-#include "itkRescaleIntensityImageFilter.h"
-
-
-int main( int argc, char *argv[] )
-{
- if( argc < 8 )
- {
- std::cerr << "Missing Parameters " << std::endl;
- std::cerr << "Usage: " << argv[0];
- std::cerr << " inputImage xi xf yi yf zi zf"<< std::endl;
- return 1;
- }
-
-
- typedef float PixelTypeFloat;
- const unsigned int Dimension = 3;
- 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();
-
- PixelTypeFloat pixelValue;
- ImageTypeFloat::IndexType pixelIndex;
-
- ImageTypeFloat::RegionType region;
- region = image->GetLargestPossibleRegion ();
-
- std::cout<<"GetImageDimensionX() :" << region.GetSize(0) << "\n" ;
- std::cout<<"GetImageDimensionY() :" << region.GetSize(1) << "\n" ;
- std::cout<<"GetImageDimensionZ() :" << region.GetSize(2) << "\n" ;
-
- std::cout<<"\nStart loops...\n" ;
- for (int k = atoi(argv[6]);k<atoi(argv[7]);k++){
- for (int i = atoi(argv[2]);i<atoi(argv[3]);i++){
- for (int j = atoi(argv[4]);j<atoi(argv[5]);j++){
- pixelIndex[0] = i; // x position
- pixelIndex[1] = j; // y position
- pixelIndex[2] = k; // y position
- pixelValue = image->GetPixel( pixelIndex );
- std::cout<< " " << pixelValue << " " ;
- }
- std::cout<< "\n" ;
- }
- std::cout<< "\n\n\n\n" ;
- }
-
- return 0;
-}
\ No newline at end of file
diff --git a/contrib/arc/conge.dat b/contrib/arc/conge.dat
deleted file mode 100644
index dc3ffaf1fa852d6861d7fc7a39ea69661388cd55..0000000000000000000000000000000000000000
--- a/contrib/arc/conge.dat
+++ /dev/null
@@ -1,13 +0,0 @@
-MeshFile conge.msh
-ElasticDomain 7 200.e9 0.3
-EdgeDisplacement 8 0 0
-EdgeDisplacement 8 1 0
-EdgeDisplacement 8 2 0
-#NodeDisplacement 20 0 0
-#NodeDisplacement 20 1 0
-#NodeDisplacement 21 1 0
-#NodeForce 21 100e6 0 0
-EdgeForce 9 100e6 0 0
-#FaceForce 7 0 -1.0e10 0
-#NodeForce 20 -1e8 0 0
-
diff --git a/contrib/arc/conge.geo b/contrib/arc/conge.geo
deleted file mode 100644
index 89cc738338b0608312660dd07b95f320059b8d83..0000000000000000000000000000000000000000
--- a/contrib/arc/conge.geo
+++ /dev/null
@@ -1,89 +0,0 @@
-unit = 1.0e-02 ;
-
-e1 = 4.5 * unit ;
-e2 = 6.0 * unit / 2.0 ;
-e3 = 5.0 * unit / 2.0 ;
-h1 = 5.0 * unit ;
-h2 = 10.0 * unit ;
-h3 = 5.0 * unit ;
-//h4 = 1.0 * unit ;
-h4 = 2.0 * unit ;
-h5 = 4.5 * unit ;
-R1 = 1.0 * unit ;
-R2 = 1.5 * unit ;
-//R2 = 1.5 * unit ;
-r = 2 * unit ;
-ccos = (-h5*R1+e2* (h5*h5+e2*e2-R1*R1)^0.5) / (h5*h5+e2*e2) ;
-ssin = ( 1.0 - ccos*ccos )^0.5 ;
-
-eps = 0.01 * unit;
-
-Lc1 = 0.001 ;
-Lc2 = 0.001 ;
-
-Point(1) = { -e1-e2, 0.0 , 0.0 , Lc1};
-Point(2) = { -e1-e2, h1 , 0.0 , Lc1};
-Point(3) = { -e3-r , h1 , 0.0 , Lc2};
-Point(4) = { -e3-r , h1+r , 0.0 , Lc2};
-Point(5) = { -e3 , h1+r , 0.0 , Lc2};
-Point(6) = { -e3 , h1+h2-eps, 0.0 , Lc1};
-Point(7) = { e3 , h1+h2, 0.0 , Lc1};
-Point(8) = { e3 , h1+r , 0.0 , Lc2};
-Point(9) = { e3+r , h1+r , 0.0 , Lc2};
-Point(10)= { e3+r , h1 , 0.0 , Lc2};
-Point(11)= { e1+e2, h1 , 0.0 , Lc1};
-Point(12)= { e1+e2, 0.0 , 0.0 , Lc1};
-Point(13)= { e2 , 0.0 , 0.0 , Lc1};
-
-Point(14)= { R1 / ssin , h5+R1*ccos, 0.0 , Lc2};
-Point(15)= { 0.0 , h5 , 0.0 , Lc2};
-Point(16)= { -R1 / ssin , h5+R1*ccos, 0.0 , Lc2};
-Point(17)= { -e2 , 0.0 , 0.0 , Lc1};
-
-Point(18)= { -R2 , h1+h3 , 0.0 , Lc2};
-Point(19)= { -R2 , h1+h3+h4, 0.0 , Lc2};
-Point(20)= { 0.0 , h1+h3+h4, 0.0 , Lc2};
-Point(21)= { R2 , h1+h3+h4, 0.0 , Lc2};
-Point(22)= { R2 , h1+h3 , 0.0 , Lc2};
-Point(23)= { 0.0 , h1+h3 , 0.0 , Lc2};
-
-Point(24)= { 0 , h1+h3+h4+R2, 0.0 , Lc2};
-Point(25)= { 0 , h1+h3-R2, 0.0 , Lc2};
-
-Line(1) = {1 ,17}; /* ux=uy=0 */
-Line(2) = {17,16};
-Circle(3) = {14,15,16};
-Line(4) = {14,13};
-Line(5) = {13,12}; /* ux=uy=0 */
-Line(6) = {12,11};
-Line(7) = {11,10};
-Circle(8) = { 8, 9,10};
-Line(9) = { 8, 7};
-Line(10) = { 7, 6}; /* T=10000 N */
-Line(11) = { 6, 5};
-Circle(12) = { 3, 4, 5};
-Line(13) = { 3, 2};
-Line(14) = { 2, 1};
-
-Line(15) = {18,19};
-Circle(16) = {21,20,24};
-Circle(17) = {24,20,19};
-Circle(18) = {18,23,25};
-Circle(19) = {25,23,22};
-Line(20) = {21,22};
-
-Line Loop(21) = {17,-15,18,19,-20,16};
-//Plane Surface(22) = {21};
-Line Loop(23) = {11,-12,13,14,1,2,-3,4,5,6,7,-8,9,10};
-Plane Surface(24) = {23,21};
-
-//Physical Line(25) = {9,1,2,3,4,5,6,7,8,11,12,13,14,15,16,17,18,19,20,10};
-//Physical Surface(26) = {22,24};
-Physical Line(9) = {11};
-Physical Line(10) = {10};
-Physical Line(8) = {1, 5};
-Physical Surface(7) = {24};
-Physical Point(20) = {2};
-Physical Point(21) = {11};
-// Physical Line(25) = {11, 12, 13, 14};
-
diff --git a/contrib/arc/conge3D.dat b/contrib/arc/conge3D.dat
deleted file mode 100644
index 412ae699d90bfd80132f57535cadbb427f8b534b..0000000000000000000000000000000000000000
--- a/contrib/arc/conge3D.dat
+++ /dev/null
@@ -1,7 +0,0 @@
-MeshFile conge3D.msh
-ElasticDomain 7 200.e9 0.3
-FaceDisplacement 8 0 0
-FaceDisplacement 8 1 0
-FaceDisplacement 8 2 0
-FaceForce 9 100e6 0 0
-VolumeForce 7 0 -1.0e10 0
diff --git a/contrib/arc/conge3D.geo b/contrib/arc/conge3D.geo
deleted file mode 100644
index 4769bb030b1e121477c88d3cea1e7bb6b80bfb92..0000000000000000000000000000000000000000
--- a/contrib/arc/conge3D.geo
+++ /dev/null
@@ -1,88 +0,0 @@
-unit = 1.0e-02 ;
-
-e1 = 4.5 * unit ;
-e2 = 6.0 * unit / 2.0 ;
-e3 = 5.0 * unit / 2.0 ;
-h1 = 5.0 * unit ;
-h2 = 10.0 * unit ;
-h3 = 5.0 * unit ;
-//h4 = 1.0 * unit ;
-h4 = 2.0 * unit ;
-h5 = 4.5 * unit ;
-R1 = 1.0 * unit ;
-R2 = 1.5 * unit ;
-//R2 = 1.5 * unit ;
-r = 2 * unit ;
-ccos = (-h5*R1+e2* (h5*h5+e2*e2-R1*R1)^0.5) / (h5*h5+e2*e2) ;
-ssin = ( 1.0 - ccos*ccos )^0.5 ;
-
-eps = 0.01 * unit;
-
-Lc1 = 0.0025 ;
-Lc2 = 0.0025 ;
-
-Point(1) = { -e1-e2, 0.0 , 0.0 , Lc1};
-Point(2) = { -e1-e2, h1 , 0.0 , Lc1};
-Point(3) = { -e3-r , h1 , 0.0 , Lc2};
-Point(4) = { -e3-r , h1+r , 0.0 , Lc2};
-Point(5) = { -e3 , h1+r , 0.0 , Lc2};
-Point(6) = { -e3 , h1+h2-eps, 0.0 , Lc1};
-Point(7) = { e3 , h1+h2, 0.0 , Lc1};
-Point(8) = { e3 , h1+r , 0.0 , Lc2};
-Point(9) = { e3+r , h1+r , 0.0 , Lc2};
-Point(10)= { e3+r , h1 , 0.0 , Lc2};
-Point(11)= { e1+e2, h1 , 0.0 , Lc1};
-Point(12)= { e1+e2, 0.0 , 0.0 , Lc1};
-Point(13)= { e2 , 0.0 , 0.0 , Lc1};
-
-Point(14)= { R1 / ssin , h5+R1*ccos, 0.0 , Lc2};
-Point(15)= { 0.0 , h5 , 0.0 , Lc2};
-Point(16)= { -R1 / ssin , h5+R1*ccos, 0.0 , Lc2};
-Point(17)= { -e2 , 0.0 , 0.0 , Lc1};
-
-Point(18)= { -R2 , h1+h3 , 0.0 , Lc2};
-Point(19)= { -R2 , h1+h3+h4, 0.0 , Lc2};
-Point(20)= { 0.0 , h1+h3+h4, 0.0 , Lc2};
-Point(21)= { R2 , h1+h3+h4, 0.0 , Lc2};
-Point(22)= { R2 , h1+h3 , 0.0 , Lc2};
-Point(23)= { 0.0 , h1+h3 , 0.0 , Lc2};
-
-Point(24)= { 0 , h1+h3+h4+R2, 0.0 , Lc2};
-Point(25)= { 0 , h1+h3-R2, 0.0 , Lc2};
-
-Line(1) = {1 ,17}; /* ux=uy=0 */
-Line(2) = {17,16};
-Circle(3) = {14,15,16};
-Line(4) = {14,13};
-Line(5) = {13,12}; /* ux=uy=0 */
-Line(6) = {12,11};
-Line(7) = {11,10};
-Circle(8) = { 8, 9,10};
-Line(9) = { 8, 7};
-Line(10) = { 7, 6}; /* T=10000 N */
-Line(11) = { 6, 5};
-Circle(12) = { 3, 4, 5};
-Line(13) = { 3, 2};
-Line(14) = { 2, 1};
-
-Line(15) = {18,19};
-Circle(16) = {21,20,24};
-Circle(17) = {24,20,19};
-Circle(18) = {18,23,25};
-Circle(19) = {25,23,22};
-Line(20) = {21,22};
-
-Line Loop(21) = {17,-15,18,19,-20,16};
-//Plane Surface(22) = {21};
-Line Loop(23) = {11,-12,13,14,1,2,-3,4,5,6,7,-8,9,10};
-Plane Surface(24) = {23,21};
-Extrude {0, 0, 0.01} {
- Surface{24};
-}
-
-
-Physical Volume(7) = {1};
-Physical Surface(10) = {101};
-Physical Surface(9) = {49};
-Physical Surface(8) = {65, 81};
-
diff --git a/contrib/arc/highlevel.h b/contrib/arc/highlevel.h
deleted file mode 100644
index e654fb24ba0362729fe00205eb072d4511480764..0000000000000000000000000000000000000000
--- a/contrib/arc/highlevel.h
+++ /dev/null
@@ -1,319 +0,0 @@
-//
-// C++ Interface: highlevel
-//
-// Description:
-//
-//
-// Author: <Eric Bechet>, (C) 2009
-//
-// Copyright: See COPYING file that comes with this distribution
-//
-//
-#ifndef highlevel_H
-#define highlevel_H
-
-#include <vector>
-#include <complex>
-#include <iostream>
-//#include <tr1/memory>
-
-#include "dofManager.h"
-#include "SVector3.h"
-#include "MElement.h"
-#include "MVertex.h"
-#include "functionSpace.h"
-
-
-
-typedef SVector3 Vector;
-class Tensor2{};
-class Tensor4{};
-
-/*
-template<class T> struct TensorialTraits
-{
- typedef T ValType;
- typedef T GradType[3];
- typedef T HessType[3][3];
- static const int nb_basis_vectors=1; // par défaut, considéré comme un scalaire
-};
-
-template<> struct TensorialTraits<double>
-{
- typedef double ValType;
- typedef Vector GradType;
- typedef Tensor2 HessType;
- static const int nb_basis_vectors=1; // scalaire
-};
-
-
-// template<> struct TensorialTraits<Vector>
-// {
-// typedef Vector ValType;
-// typedef Tensor2 GradType;
-// typedef Tensor3 HessType;
-// static const int nb_basis_vectors=3; // trois vecteurs de base linéairement indépendants.
-// };
-*/
-
-template<class T> class Function // Fonction au sens EF du terme.
- // renvoie valeur , gradient, hessien, etc...pour un element donné et un point donné
-{
-public:
-// typedef std::tr1::shared_ptr<Function<T> > FunctionPtr;
- virtual ~Function(){}
- virtual void GetVal (double uvw[3],MElement *e, T& Val)const=0;
- virtual void GetGrad(double uvw[3],MElement *e,typename TensorialTraits<T>::GradType &Grad) const =0;
- virtual void GetHess(double uvw[3],MElement *e,typename TensorialTraits<T>::HessType &Hess)const =0;
- virtual void GetVal (double uvw[][3],MElement *e, T Val[],int n) { for (int i=0;i<n;++i) GetVal(uvw[i],e,Val[i]); } // par defaut
- virtual void GetGrad(double uvw[][3],MElement *e,typename TensorialTraits<T>::GradType Grad[],int n) const { for (int i=0;i<n;++i) GetGrad(uvw[i],e,Grad[i]); }
- virtual void GetHess(double uvw[][3],MElement *e,typename TensorialTraits<T>::HessType Hess[],int n)const { for (int i=0;i<n;++i) GetHess(uvw[i],e,Hess[i]); }
-};
-
-
-template<class T> class LagrangeShapeFunction: public Function<T>
-{
- public:
- virtual void GetVal(double uvw[],MElement *e, T& Val)
- {
-// double s[100];
-// _e->getShapeFunctions(uvw[0], uvw[1], uvw[2], s);
-
- }
-};
-
-
-
-
-
-/*
-class SpaceBase // renvoie des clefs des dofs
-{
-protected:
-int _iField;
-public:
- SpaceBase(int iField=0):_iField(iField){};
- virtual ~SpaceBase(){};
- virtual int getNumberDofs(MElement *e) const = 0 ;
- virtual void getDofs(MElement *e,const Dof *vecD) const {} //=0;
- virtual int getId(void) const {return _iField;};
-};
-
-template<class T> class Space : public SpaceBase // renvoie des clefs de dofs et des fonctions de formes
-{
-
-public:
- Space(int iField=0):SpaceBase(iField){};
- virtual ~Space(){};
- virtual void getDofsandSFs(MElement *e,const Dof *VecD, const Function<T>* vecDSF) {} //=0;
- virtual void getSFs(const Function<T>* vecDSF) {} //=0;
-};
-
-
-template<class T> class SpaceLagrange : public Space<T> // Lagrange ... 1 dof / node / basis vector
-{
-private:
-
- Dof getLocalDof(MElement *e, int i) const
- {
- int iComp = i / e->getNumVertices();
- int ithLocalVertex = i % e->getNumVertices();
- return Dof(e->getVertex(ithLocalVertex)->getNum(),
- Dof::createTypeWithTwoInts(iComp,Space<T>::getId()));
- }
-
-public:
- SpaceLagrange(int iField=0):Space<T>(iField){};
- virtual ~SpaceLagrange(){};
- virtual void getDofsandSFs(MElement *e,std::vector<std::pair< Dof, Function<T>* > > &vecDFF){}
- virtual void getSFs(std::vector<std::pair< Dof, Function<T>* > > &vecDFF)
- {
-
- };
- virtual void getDofs(MElement *e,std::vector<Dof> &vecD) const
- {
- int ndofs= e->getNumVertices()*TensorialTraits<T>::nb_basis_vectors;
- for (int i=0;i<ndofs;++i)
- vecD.push_back(getLocalDof(e,i));
- }
- virtual int getNumberDofs(MElement *e) const {return 0;}
-};
-
-
-template<class T> class SpaceXfem : public Space<T> // approximation xfem (spacebase + enrich)
-{
- int _iEnrich;
- Space<T> &_SpaceBase;
- Function<double> &enrich;
- std::unary_function<MVertex*,bool> &test;
- Dof getLocalDof(MElement *e, int i) const
- {
- int iComp = i / e->getNumVertices();
- int ithLocalVertex = i % e->getNumVertices();
- MVertex * v= e->getVertex(ithLocalVertex);
- if (test(v))
- return Dof(e->getVertex(ithLocalVertex)->getNum(),
- Dof::createTypeWithTwoInts(iComp, _iEnrich));
- }
-public:
- SpaceXfem(int iEnrich,Space<T>& SpaceBase):_iEnrich(iEnrich), _SpaceBase(SpaceBase){};
- virtual void getDofsandSFs(MElement *e,std::vector<std::pair< Dof, Function<T>* > > &vecDFF){}
- virtual void getSFs(std::vector<std::pair< Dof, Function<T>* > > &vecDFF){};
- virtual void getDofs(MElement *e,std::vector<Dof> &vecD) const
- {
- _SpaceBase.getDofs(e,vecD);
- }
-};
-
-*/
-
-/*
-template<class T> class Field : public Function<T> // renvoie des valeurs de champ (ff*valeurs dofs), gradient , etc...
-{
-public:
- Field(){};
- virtual ~Field(){};
-};
-*/
-
-
-class FormBilinearBase
-{
- public:
- template <class T> void allocate(T *p) {p=new T[10];std::cout << "10" << std::endl;}
- template <class T> void get(T* tab, T *p, int ndx) {p=&tab[ndx];}
- virtual void getDofs(MElement *e)=0;
- virtual void getFuncs(MElement *e,const double uvw[3])=0;
- virtual void getGradFuncs(MElement *e,const double uvw[3])=0;
-};
-
-template <> void FormBilinearBase::allocate(void *p) {p=NULL; std::cout << "null" << std::endl;}
-template <> void FormBilinearBase::get(void *tab,void *p, int ndx) {p=NULL;}
-
-template <class Term,class SpaceL,class SpaceR> class FormBilinear : public FormBilinearBase
- // Renvoie des matrices élémentaires (ff)
- // Accessoirement stocke des pointeurs vers les termes pour chaque element
- // Doit etre initialisée AVANT toute opération (pour l'allocation)
- // en principe ce truc ne devrait pas ^etre reimplemente
- // il devrait donc dependre d'un parametre template TermBilinear
- // elle sait "integrer" dans un elemeent
-{
-
- typedef typename Term::Stype Stype; // le truc a stocker
- typedef typename Term::Datamat Datamat; // type elementaire
- Term instance; // une instance du terme. E.g. utile pour toute initialisation (calcul tenseur materiel, etc...)
- Stype *p;
- MElement *e;
- Function<double> *a;
- Function<double> *b;
- double uvw[3];
- fullMatrix<Datamat> mat;
- SpaceL *_spacel;
- SpaceR *_spacer;
- std::vector<Dof> Dofsl;
- std::vector<Dof> Dofsr;
- std::vector<Datamat> Funcl;
- std::vector<Datamat> Funcr;
-
- public:
- FormBilinear(SpaceL &sl,SpaceR &sr) {}
- void getMatrix(fullMatrix<Datamat> *v) {v=&mat;}
- void func(void) { allocate(p);}
- void func2(void) { }
- void getDofs(MElement *e)
- {
- _spacel->getKeys(e,Dofsl);
- _spacer->getKeys(e,Dofsr);
- };
- void getFuncs(MElement *e,const double uvw[3]) {};
- void getGradFuncs(MElement *e,const double uvw[3]) {};
- void Init(MElement *e) { getDofs(e);} // called once for each element
- void Accumulate(MElement *e,const double uvw[3]) // called for every GP
- {
-// for (int i)
-// for (int j)
-// mat(i,j)+=instance.getTerm(uvw,*e,*a,*b,NULL,this);
- instance.Init(e,this); // called once for each GP
- if (instance.NeedFuncs()) getFuncs(e,uvw);
- if (instance.NeedGradFuncs()) getGradFuncs(e,uvw);
- Datamat result;
- instance.getTerm(uvw,e,NULL,*a,*b,result); // called for every SF combination (times every GP, times every element)
- }
-};
-
-//template <class Term,class SpaceL,class datamat=double> class FormBilinear<Term,SpaceL,SpaceR> : public FormBilinearBase {}; cas symétrique
-
-
-class FormLinear{}; // renvoie des vecteurs élémentaires
- // on devrait pouvoir construire une forme lin à partir d'une forme bilin pour les pb "matrix free"
- // idem FormBilin
-
-class FormZero{}; // renvoie des scalaires ex. resultat d'une integration
-
-//algorithmes
-/*void Construct(FormBilinear &B,Region &R,Integrator &I);
-void Assemble(FormBilinear &B,Region &R,Assembler &A, Integrator &I);
-void Construct(FormLinear &L,Region &R,Integrator &I);
-void Assemble(FormLinear &L,Region &R,Assembler &A,Integrator &I);
-void Construct(FormZero &Z,Region &R,Integrator &I);
-void Assemble(FormZero &Z,Region &R,Assembler &A,Integrator &I);
-*/
-
-class TermBilinearBase
-{
-public:
- virtual void GetTerm(MElement *e) {}
-};
-
-
-template<class datamat, class T1,class T2> class TermBilinear : public TermBilinearBase // terme associe a un "element" / pt de gauss (contribution élémentaire)
- // typiquement celui ci stoque ce qui doit etre stocke. C'est la base configurable du code
- // on doit associer cela a un allocateur qui renvoie un pointeur sur ce truc
- // Soit tous les elements/pts de gauss ont le meme terme , allocateur unique (pas de stockage aux pts de
- // gauss par exemple
- // soit ils ont qqc a stockuer et sont tous distincts
- // on peut utiliser des membres statiques pour ce qui est constant pour tous les instances
-{
-public:
- typedef void Stype;
- typedef datamat Datamat;
-// TermBilinear(){};
- virtual void Init(MElement *e, FormBilinearBase * caller) //called once for every element
- {
- }
- virtual bool NeedGradFuncs(void) { return false ;}// Query if getTerm will use the gradient of the SF
- virtual bool NeedFuncs(void) { return false ;}// Query if it will need the value of the SF
-
- virtual void getTerm(const double uvw[],MElement *e,Stype* info, Function<T1> &SF,Function<T2> &TF, datamat &res ) // called for every gauss point
- {
- };
- virtual void Update(const double uvw[],MElement &e,Function<T1> &SF,Function<T2> &TF) {};
-// toute fonctios utiles . prevoir un algorithme ad hoc pour l'appliquer sur ts les pts de gauss ...
-};
-
-
-class TermBilinearMeca : public TermBilinear<double, double,double>
-{
-public:
- virtual bool NeedGradFuncs(void) { return true ;}// Query if getTerm will use the gradient of the SF
-// typedef Tensor Stype;
-// TermBilinear(){};
-// virtual double getTerm(double uvw[],MElement &e,Function<double> &SF,Function<double> &TF,Stype* info) {};
-// virtual void Update(double uvw[],MElement &e,Function<T1> &SF,Function<T2> &TF) {};
-// toute fonctios utiles . prevoir un algorithme
-};
-
-
-class TermBilinearMecaNL : public TermBilinearMeca
-{
-public:
- typedef Tensor2 Stype;
-// TermBilinear(){};
-// virtual double getTerm(double uvw[],MElement &e,Function<double> &SF,Function<double> &TF,Stype* info) {};
-// virtual void Update(double uvw[],MElement &e,Function<T1> &SF,Function<T2> &TF) {};
-// toute fonctios utiles . prevoir un algorithme
-};
-
-
-
-#endif // highlevel_H
diff --git a/contrib/arc/mainElasticity.cpp b/contrib/arc/mainElasticity.cpp
deleted file mode 100644
index 24194ca1b15681badf54483f2a036c30f3e1f7b0..0000000000000000000000000000000000000000
--- a/contrib/arc/mainElasticity.cpp
+++ /dev/null
@@ -1,120 +0,0 @@
-#include "Gmsh.h"
-#include "elasticitySolver.h"
-#include "PView.h"
-#include "PViewData.h"
-#include "highlevel.h"
-#include "groupOfElements.h"
-#include <iterator>
-#include "function.h"
-#include "fullMatrix.h"
-
-class functionAdd : public function
-{
- private:
- std::vector<std::string> strvec;
- class data : public dataCacheDouble
- {
- private:
- std::vector<dataCacheDouble *> dcvec;
- public:
- data(const functionAdd * fm,dataCacheMap *m) :
- dataCacheDouble(*m,m->getNbEvaluationPoints(),1)
- {
- for (int i=0;i<fm->strvec.size();++i)
- dcvec.push_back(&(m->get(fm->strvec[i],this)));
- }
- void _eval()
- {
- _value(0, 0)=0;
- for (int i=0;i<dcvec.size();++i) { _value(0, 0) += (*(dcvec[i]))(0,0);}
- }
- ~data()
- {
- }
- };
- public:
- dataCacheDouble *newDataCache(dataCacheMap *m)
- {
- return new data(this,m);
- }
- functionAdd(std::string _a,std::string _b) {strvec.push_back(_a);strvec.push_back(_b);}
- void addNewTerm(std::string _a) { strvec.push_back(_a);}
-};
-
-
-
-
-int main (int argc, char* argv[])
-{
-
- if (argc != 2){
- printf("usage : elasticity input_file_name\n");
- return -1;
- }
-#if 0
- fullMatrix<double> a(1,1);
- a(0,0)=1.0;
- fullMatrix<double> b(1,1);
- b(0,0)=2.0;
- fullMatrix<double> c(1,1);
- c(0,0)=4.0;
- fullMatrix<double> d(1,1);
- d(0,0)=8.0;
- fullMatrix<double> e(1,1);
- e(0,0)=16.0;
-
- std::cout << argc << std::endl;
- functionConstant fc_a(&a);
- functionConstant fc_b(&b);
- functionConstant fc_c(&c);
- functionConstant fc_d(&d);
- functionConstant fc_e(&e);
-
- functionAdd fa_ab(fc_a.getName(),fc_b.getName());
- functionAdd fa_cd(fc_c.getName(),fc_d.getName());
- functionAdd fa_abcd("a+b","c+d");
-
- function::add("a+b", &fa_ab );
- function::add("c+d", &fa_cd);
- function::add("a+b+c+d", &fa_abcd);
-
- dataCacheMap m(1);
- std::cout << "start" << std::endl;
- dataCacheDouble &dc_a=m.get(fc_a.getName());
- dataCacheDouble &dc_abcd=m.get("a+b+c+d");
- fa_abcd.addNewTerm(fc_e.getName());
- std::cout << "nb depsabcd=" << dc_abcd.howManyDoIDependOn() << std::endl;
- std::cout << "nb depsa=" << dc_a.howManyDependOnMe() << std::endl;
- std::cout << "a" << std::endl;
- std::cout << dc_a(0,0) << std::endl;
- std::cout << "a+b+c+d" << std::endl;
- std::cout << dc_abcd(0,0) << std::endl;
- std::cout << "dca.set(b)" << std::endl;
- dc_a.set(b);
- std::cout << "a+b+c+d" << std::endl;
- std::cout << dc_abcd(0,0) << std::endl;
- return(0);
-#endif
-
- GmshInitialize(argc, argv);
- // globals are still present in Gmsh
-
- // instanciate a solver
- elasticitySolver 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();
-}
diff --git a/contrib/arc/mainImageSolver.cpp b/contrib/arc/mainImageSolver.cpp
deleted file mode 100644
index b5ce919b8d66c5ba4d1b909621f526d65f95cac4..0000000000000000000000000000000000000000
--- a/contrib/arc/mainImageSolver.cpp
+++ /dev/null
@@ -1,205 +0,0 @@
-//
-// 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 "GModel.h"
-#include "Gmsh.h"
-#include "elasticitySolver.h"
-#include "PView.h"
-#include "PViewData.h"
-
-#include "groupOfElements.h"
-#include <iterator>
-
-
-#include "Classes/ImageSolver.cpp"
-
-#include "Classes/QuadtreeLSImage.cpp"
-#include "Classes/OctreeLSImage.cpp"
-
-
-
-int main( int argc, char *argv[] )
-{
- if( argc < 5 )
- {
- std::cerr << "Missing Parameters " << std::endl;
- std::cerr << "Usage: " << argv[0];
- std::cerr << " inputLevelSetImage MeshSizeMax MeshSizeMin input.dat"<< std::endl;
- return 1;
- }
-
- bool simplex = 1;
- double facx = 1;
- double facy = 1;
- double facz = 1;
-
- GmshInitialize(argc, argv);
-
- typedef float PixelTypeFloat;
- const unsigned int Dimension = 3;
- typedef itk::Image< PixelTypeFloat, Dimension > ImageTypeFloat;
- typedef itk::ImageFileReader< ImageTypeFloat > ImageReaderTypeFloat;
- ImageReaderTypeFloat::Pointer reader = ImageReaderTypeFloat::New();
-
- typedef itk::Image< PixelTypeFloat, 2 > ImageTypeFloat2D;
-
-
- reader->SetFileName(argv[1]);
-
- ImageTypeFloat::Pointer image3D = reader->GetOutput();
- image3D->Update();
-
- ImageTypeFloat::RegionType region;
- region = image3D->GetLargestPossibleRegion ();
- std::cout<<"\nImage dimensions : " << region.GetSize(0) << " x " << region.GetSize(1) << " x " << region.GetSize(2)<<"\n";
-
-
- if (region.GetSize(2)==1)
- {
- typedef itk::ImageFileReader< ImageTypeFloat2D > ImageReaderTypeFloat2D;
- ImageReaderTypeFloat2D::Pointer reader2D = ImageReaderTypeFloat2D::New();
- ImageTypeFloat2D::Pointer image2D = reader2D->GetOutput();
- image2D->Update();
- QuadtreeLSImage tree(image2D);
-
-
- int sizemax = atoi(argv[2]);
- int sizemin = atoi(argv[3]);
-
- // Mesh with conditions on mesh size
- tree.Mesh(sizemax,sizemin);
-
- // Smoothing mesh to avoid too much generation between adjacent elements
- bool statut=false;
- int k = 0;
- std::cout<<"\nLeaf Number : "<< (tree.GetLeafNumber())<<"\n";
- while((!statut) & (k < 20)){
- statut = tree.Smooth();
- std::cout<<"\nk : "<< k;
- std::cout<<"\nsmoothed : " << (int)statut;
- std::cout<<"\nLeaf Number : "<< (tree.GetLeafNumber())<<"\n";
- k++;
- }
-
- // Create GModel with the octree mesh
- GModel *m = tree.CreateGModel(simplex,facx,facy,facz,sizemax, sizemin);
-
- // Write a .msh file with the mesh
- std::string ModelName = "TreeMesh.msh" ;
- m->writeMSH(ModelName,2.1,false,false);
-
- // Write a .msh file with the level set values as postview data
- PView *pv = tree.CreateLSPView(m);
- bool useadapt = true;
- pv->getData(useadapt)->writeMSH("LSPView.msh", false);
-
-
- ImageSolver mySolver(1000);
-
- std::vector< float >* ListLSValue = tree.getListLSValue();
- std::map< int ,std::vector<int> > * ListHangingNodes = tree.getListHangingNodes();
- mySolver.setMesh(m,2,*ListLSValue,*ListHangingNodes);
-
- int ok;
- std::cout<<"Verifier les conditions limites \n";
- //std::cin>>ok;
-
- mySolver.readInputFile(argv[4]);
-
- //solve the problem
- mySolver.solve() ;
-
- 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;
- }
- if (region.GetSize(2)>1)
- {
-
- OctreeLSImage tree(image3D);
-
-
- int sizemax = atoi(argv[2]);
- int sizemin = atoi(argv[3]);
-
- // Mesh with conditions on mesh size
- tree.Mesh(sizemax,sizemin);
-
- // Smoothing mesh to avoid too much generation between adjacent elements
- bool statut=false;
- int k = 0;
- std::cout<<"\nLeaf Number : "<< (tree.GetLeafNumber())<<"\n";
- while((!statut) & (k < 20)){
- statut = tree.Smooth();
- std::cout<<"\nk : "<< k;
- std::cout<<"\nsmoothed : " << (int)statut;
- std::cout<<"\nLeaf Number : "<< (tree.GetLeafNumber())<<"\n";
- k++;
- }
-
- // Create GModel with the octree mesh
- GModel *m = tree.CreateGModel(simplex,facx,facy,facz,sizemax, sizemin);
-
- // Write a .msh file with the mesh
- std::string ModelName = "TreeMesh.msh" ;
- m->writeMSH(ModelName,2.1,false,false);
-
- PView *pv;
- // Write a .msh file with the level set values as postview data
- //pv = tree.CreateLSPView(m);
-// bool useadapt = true;
-// pv->getData(useadapt)->writeMSH("LSPView.msh", false);
-
-
-
- ImageSolver mySolver(1000);
-
- std::vector< float >* ListLSValue = tree.getListLSValue();
- std::map< int ,std::vector<int> > * ListHangingNodes = tree.getListHangingNodes();
- mySolver.setMesh(m,3,*ListLSValue,*ListHangingNodes);
-
- int ok;
- std::cout<<"Verifier les conditions limites \n";
- //std::cin>>ok;
-
- mySolver.readInputFile(argv[4]);
-
- // solve the problem
- mySolver.solve() ;
-//
- 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();
-
- std::cout<<"\n";
-
- return 0;
-}