From 6a865af952457f94ded7d6870ebb795b1f6f57b4 Mon Sep 17 00:00:00 2001
From: Boris Sedji <sedji.boris@hotmail.com>
Date: Tue, 2 Mar 2010 12:40:05 +0000
Subject: [PATCH]
---
contrib/arc/Classes/FuncGradDisc.h | 6 +-
contrib/arc/Classes/ImageSolver.cpp | 377 +++++
contrib/arc/Classes/ImageSolver.h | 54 +
contrib/arc/Classes/OctreeLSImage.cpp | 2036 +++++++++++++----------
contrib/arc/Classes/OctreeLSImage.h | 6 +
contrib/arc/Classes/OctreeSolver.cpp | 4 +-
contrib/arc/Classes/QuadtreeLSImage.cpp | 1344 ++++++++-------
contrib/arc/Classes/QuadtreeLSImage.h | 37 +-
contrib/arc/XFEMInclusion.cpp | 114 ++
contrib/arc/mainImageSolver.cpp | 205 +++
10 files changed, 2620 insertions(+), 1563 deletions(-)
create mode 100644 contrib/arc/Classes/ImageSolver.cpp
create mode 100644 contrib/arc/Classes/ImageSolver.h
create mode 100644 contrib/arc/XFEMInclusion.cpp
create mode 100644 contrib/arc/mainImageSolver.cpp
diff --git a/contrib/arc/Classes/FuncGradDisc.h b/contrib/arc/Classes/FuncGradDisc.h
index 375497f6a0..dc09e62df6 100644
--- a/contrib/arc/Classes/FuncGradDisc.h
+++ b/contrib/arc/Classes/FuncGradDisc.h
@@ -63,13 +63,13 @@ class FuncGradDisc : public simpleFunction<double> {
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";
+ //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";
+ //std::cout<<"val f :" << f << "\n";
return f;
diff --git a/contrib/arc/Classes/ImageSolver.cpp b/contrib/arc/Classes/ImageSolver.cpp
new file mode 100644
index 0000000000..8bb0a81019
--- /dev/null
+++ b/contrib/arc/Classes/ImageSolver.cpp
@@ -0,0 +1,377 @@
+//
+// 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
new file mode 100644
index 0000000000..c1de3ab7ee
--- /dev/null
+++ b/contrib/arc/Classes/ImageSolver.h
@@ -0,0 +1,54 @@
+//
+// 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/OctreeLSImage.cpp b/contrib/arc/Classes/OctreeLSImage.cpp
index f6839e02ae..e248d77187 100644
--- a/contrib/arc/Classes/OctreeLSImage.cpp
+++ b/contrib/arc/Classes/OctreeLSImage.cpp
@@ -14,179 +14,174 @@
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];
- }
+ 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;
+ }
- 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];
- }
+ 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;
+ // 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";
+ std::cout<<"size root :"<<_root->BoxSize()<<"\n";
}
void OctreeLSImage::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);
- };
+ 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;
+ _data = data;
+ _parent = parent;
+ _children[0] = NULL;
+ _octree = octree;
}
bool Box::HaveChildren(){
- if (this->_children[0]==NULL) return false;
- else return true;
+ 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 ();
+ 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++){
+ // 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[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[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
+ 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;
+ pixelValue = image->GetPixel( pixelIndex );
+ data->LevelSetValue[i][j][k] = pixelValue;
- }
- }
- }
+ }
+ }
+ }
}
int Box::BoxSize(){
- int SizeX;
- int SizeY;
- int SizeZ;
- int size;
+ 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];
+ 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 (SizeX<SizeY) size = SizeX;
+ else size = SizeY;
- if (SizeZ<size) size = SizeZ;
+ if (SizeZ<size) size = SizeZ;
- return size;
+ 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;
+ 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;
}
@@ -195,29 +190,34 @@ bool Box::IsHomogeneous(){
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
+ }
- 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();
+ // else we divide...
+ this->Divide();
- for (int n = 0;n<8;n++){
- _children[n]->Mesh(maxsize,minsize);
- }
+ for (int n = 0;n<8;n++){
+ _children[n]->Mesh(maxsize,minsize);
+ }
}
-void Box::PrintLevelSetValue(){
+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];
- }
- }
- }
+ 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];
+ }
+ }
+ }
}
@@ -225,570 +225,806 @@ void Box::PrintLevelSetValue(){
GModel *OctreeLSImage::CreateGModel(bool simplex, double facx, double facy, double facz,int & sizemax,int & sizemin){
- this->FillMeshInfo(sizemin);
-
- 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++){
-// 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 | _ListNodes[(*ite)[j]-1][2] > zlimit ){
- inImage = false;
+ 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++;
}
- }
-
- 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();
- }
+ }
+ 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
}
+ }
- 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)[3]);
- indices.push_back((*ite)[6]);
- 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++;
+ 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);
- // --- 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]]->z() > (zlimit-eps)*facz)
- {
- front_6_indices.push_back(indices[j]);
- }
- if (vertexMap[indices[j]]->x() > (xlimit-eps)*facx)
- {
- front_4_indices.push_back(indices[j]);
- }
- if (vertexMap[indices[j]]->y() < eps*facy)
- {
- front_2_indices.push_back(indices[j]);
- }
- }
+ std::cout<<"getNumScalars :"<< pv->getData(useadapt)->getNumScalars()<<"\n";
- 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_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_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++;
- }
+//
+// 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;
- // --- Domain element tetraedre 2 ---
- indices.clear();
- indices.push_back((*ite)[6]);
- indices.push_back((*ite)[3]);
- indices.push_back((*ite)[1]);
- 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++;
+}
+PView *OctreeLSImage::CreateLSPView(GModel* m){
- // --- Domain element tetraedre 3 ---
- indices.clear();
- indices.push_back((*ite)[3]);
- indices.push_back((*ite)[4]);
- indices.push_back((*ite)[0]);
- indices.push_back((*ite)[1]);
- 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);
+ 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;
- // --- 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]]->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 (vertexMap[indices[j]]->z() < eps*facz)
- {
- front_3_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_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++;
- }
- // --- Domain element tetraedre 4 ---
- indices.clear();
- indices.push_back((*ite)[6]);
- indices.push_back((*ite)[3]);
- indices.push_back((*ite)[2]);
- indices.push_back((*ite)[1]);
- 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++;
+void OctreeLSImage::FillHangingNodes(){
- // --- 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() < eps*facx)
- {
- front_1_indices.push_back(indices[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_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++;
- }
+ _root->FillHangingNodes(_ListNodes,_ListElements,_HangingNodes);
- 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++;
- }
+void Box::FillHangingNodes(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements,std::map< int ,std::vector<int> > &HangingNodes){
- // --- Domain element tetraedre 5 ---
- indices.clear();
- indices.push_back((*ite)[4]);
- indices.push_back((*ite)[1]);
- indices.push_back((*ite)[6]);
- indices.push_back((*ite)[5]);
- 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++;
+ std::vector<std::vector<int> >::iterator it;
- // --- 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++)
+ 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)
{
- 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 (vertexMap[indices[j]]->z() < eps*facz)
- {
- front_3_indices.push_back(indices[j]);
- }
+ 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 (front_4_indices.size() == 3)
+ if (Leafs.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++;
+ 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 (front_5_indices.size() == 3)
+ if (Leafs.size()==6)
{
- 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++;
+ 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 (front_3_indices.size() == 3)
+ if (Leafs.size()==7)
{
- 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++;
+ 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;
}
- }
- ite++; // next octree leaf
+ it++;
+ k++;
}
- std::cout<<"k :"<<k<<"\n";
- }
-
- 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::FillMeshInfo(int & pas){
- if (!_ListNodes.empty()){
- return;
- }
+ if (!_ListNodes.empty()){
+ return;
+ }
- _root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue,pas);
+ _root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue,pas);
}
@@ -800,79 +1036,85 @@ void OctreeLSImage::FillMeshInfo(int & pas){
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::cout<<"\nxyzpas :"<< xpas << " " << ypas << " " << zpas;
-
- 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();
- _octree->GetRoot()->FillElementsNode(ListElements);
+ 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;
+ // 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];
@@ -888,21 +1130,21 @@ void Box::FillElementsNode(std::vector< std::vector<int> > &ListElements){
// ElementNodes[6] = ElementNodes[7];
// ElementNodes[7] = temp;
- ListElements.push_back(ElementNodes);
+ ListElements.push_back(ElementNodes);
// std::cout<<"\n Nodes :";
// for (int i = 0 ; i < 8 ; i++)
// std::cout<<" " << ElementNodes[i];
- return;
- }
+ return;
+ }
- for (int n = 0 ; n < 8 ; n++){
- this->_children[n]->FillElementsNode(ListElements);
- }
+ for (int n = 0 ; n < 8 ; 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){
+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()) {
@@ -1008,19 +1250,19 @@ void Box::GetLeafElements(std::vector< std::vector<int> > &ListElements){
void OctreeLSImage::SetLeafNumber(){
_LeafNumber = 0;
- _root->CountLeafNumber(_LeafNumber);
+ _root->CountLeafNumber(_LeafNumber);
}
void Box::CountLeafNumber(int &LeafNumber){
- if (this->GetChild(0)==NULL) {
- LeafNumber = LeafNumber + 1;
- return;
- }
+ if (this->GetChild(0)==NULL) {
+ LeafNumber = LeafNumber + 1;
+ return;
+ }
- for (int n = 0 ; n < 8 ; n++){
- this->_children[n]->CountLeafNumber(LeafNumber);
- }
+ for (int n = 0 ; n < 8 ; n++){
+ this->_children[n]->CountLeafNumber(LeafNumber);
+ }
}
@@ -1030,132 +1272,132 @@ void Box::Divide(){
// box divide in 8 boxes whit limite xl = (xmin + xmax)/2 ; yl = (ymin + ymax)/2 ; zl = (zmin + zmax)/2
- int n;
+ int n;
- // Box 1 : < xl ; < yl ; < zl
+ // Box 1 : < xl ; < yl ; < zl
- n = 0;
- BoxData* data1 = new BoxData;
+ 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;
+ 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());
+ 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
+ // Box 2 : > xl ; < yl ; < zl
- n = 1;
- BoxData* data2 = new BoxData;
+ 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;
+ 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());
+ 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
+ // Box 3 : < xl ; > yl ; < zl
- n = 2;
- BoxData* data3 = new BoxData;
+ 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;
+ 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());
+ 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
+ // Box 4 : > xl ; > yl ; < zl
- n = 3;
- BoxData* data4 = new BoxData;
+ 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;
+ 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());
+ 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
+ // Box 5 : < xl ; < yl ; > zl
- n = 4;
- BoxData* data5 = new BoxData;
+ 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];
+ 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());
+ 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
+ // Box 6 : > xl ; < yl ; > zl
- n = 5;
- BoxData* data6 = new BoxData;
+ 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];
+ 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());
+ 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
+ // Box 7 : < xl ; > yl ; > zl
- n = 6;
- BoxData* data7 = new BoxData;
+ 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];
+ 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());
+ 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
+ // Box 8 : > xl ; > yl ; > zl
- n = 7;
- BoxData* data8 = new BoxData;
+ 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];
+ 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());
+ 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();
+ return _root->Smooth();
}
@@ -1163,86 +1405,86 @@ bool OctreeLSImage::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;
- }
+ 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;
+ 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()){
+ if (!this->HaveChildren()){
Box* parent = this;
BoxData* data;
- while(parent->GetParent()!=NULL){
- data = parent->GetParent()->GetData();
- parent=parent->GetParent();
+ while (parent->GetParent()!=NULL){
+ data = parent->GetParent()->GetData();
+ parent=parent->GetParent();
}
- Leafs.push_back(this);
- return;
- }
+ 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[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[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[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[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[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[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[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);
+ 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
index 3b732115f6..f083594284 100644
--- a/contrib/arc/Classes/OctreeLSImage.h
+++ b/contrib/arc/Classes/OctreeLSImage.h
@@ -49,11 +49,14 @@ class OctreeLSImage
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 :
@@ -74,6 +77,8 @@ class OctreeLSImage
PView* CreateLSPView(GModel* m);
int* GetSize(){return _ImageSize;}
+ std::vector< float >* getListLSValue(){return &_ListLSValue;}
+ std::map< int ,std::vector<int> >* getListHangingNodes(){return &_HangingNodes;}
};
@@ -117,6 +122,7 @@ class Box{
// 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
index 2401e6f013..f34998b3b2 100644
--- a/contrib/arc/Classes/OctreeSolver.cpp
+++ b/contrib/arc/Classes/OctreeSolver.cpp
@@ -8,8 +8,6 @@
//
//
-
-
#include <string.h>
#include "GmshConfig.h"
#include "elasticitySolver.h"
@@ -248,7 +246,7 @@ void OctreeSolver::solve(){
//_TagEnrichedVertex.clear();
_EnrichComp.insert(0);
_EnrichComp.insert(1);
- _EnrichComp.insert(2);
+ //_EnrichComp.insert(2);
_funcEnrichment = new FuncGradDisc(_ls,pModel);
diff --git a/contrib/arc/Classes/QuadtreeLSImage.cpp b/contrib/arc/Classes/QuadtreeLSImage.cpp
index 9a2c0e457e..4a49995bb2 100644
--- a/contrib/arc/Classes/QuadtreeLSImage.cpp
+++ b/contrib/arc/Classes/QuadtreeLSImage.cpp
@@ -14,590 +14,645 @@
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 ();
+ 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
- // Size will be a power of two larger than the image size
+ int size[3];
+ size[0] = 1;
+ size[1] = 1;
- int size[3];
- size[0] = 1;
- size[1] = 1;
+ _ImageSize[0] = region.GetSize(0);
+ _ImageSize[1] = region.GetSize(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;
- }
-
- BoxData *data = new BoxData;
+ while (size[0]<region.GetSize(0)){
+ size[0] = size[0]*2;
+ }
+ while (size[1]<region.GetSize(1)){
+ size[1] = size[1]*2;
+ }
- // The window of the image root
+ Box2DData *data = new Box2DData;
- 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];
+ 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 Box(data,NULL,this);
- _root->FillLevelSetValue(image,data);
- _image = image;
- _LeafNumber = 0;
+ }
+
+ // 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
+ // taillemax //nombre de pixels max dans une direction
- if (this->_root->HaveChildren()){
- return;
- }
- if(_root->BoxSize()>maxsize || !_root->IsHomogeneous()){
- _root->Mesh( maxsize, minsize);
- };
+ if (this->_root->HaveChildren()){
+ return;
+ }
+ if (_root->BoxSize()>maxsize || !_root->IsHomogeneous()){
+ _root->Mesh( maxsize, minsize);
+ };
}
-Box::Box(BoxData* data, Box* parent,QuadtreeLSImage* quadtree){
+Box2D::Box2D(Box2DData* data, Box2D* parent,QuadtreeLSImage* quadtree){
- _data = data;
- _parent = parent;
- _children[0] = NULL;
- _quadtree = quadtree;
+ _data = data;
+ _parent = parent;
+ _children[0] = NULL;
+ _quadtree = quadtree;
}
-bool Box::HaveChildren(){
+bool Box2D::HaveChildren(){
- if (this->_children[0]==NULL) return false;
- else return true;
+ if (this->_children[0]==NULL) return false;
+ else return true;
}
-void Box::FillLevelSetValue(itk::Image< float, 2 >::Pointer image, BoxData *data){
+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 ();
+ 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++){
+ // 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[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[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;
+ pixelValue = image->GetPixel( pixelIndex );
+ data->LevelSetValue[i][j] = pixelValue;
- }
- }
+ }
+ }
}
-int Box::BoxSize(){
+int Box2D::BoxSize(){
- int SizeX;
- int SizeY;
- int size;
+ int SizeX;
+ int SizeY;
+ int size;
- SizeX = this->_data->boundingbox[1] - this->_data->boundingbox[0];
- SizeY = this->_data->boundingbox[3] - this->_data->boundingbox[2];
+ 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;
+ if (SizeX<SizeY) size = SizeX;
+ else size = SizeY;
- return size;
+ return size;
}
-bool Box::IsHomogeneous(){
+bool Box2D::IsHomogeneous(){
- float OneValue = this->_data->LevelSetValue[0][0];
+ 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;
+ 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 Box::Mesh(int & maxsize, int & minsize){
+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
- }
+ 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();
+ // else we divide...
+ this->Divide();
- for (int n = 0;n<4;n++){
- _children[n]->Mesh(maxsize,minsize);
- }
+ for (int n = 0;n<4;n++){
+ _children[n]->Mesh(maxsize,minsize);
+ }
}
-void Box::PrintLevelSetValue(){
+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];
- }
- }
+ 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,int & sizemax,int & sizemin){
-
- double eps = 1e-6;
- this->FillMeshInfo(sizemin);
-
- 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++;
+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++)
- {
+ 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;
+ 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
}
- }
-
- //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++;
- }
+ std::cout<<"numElement size :"<<numElement.size()<<"\n";
+ std::cout<<"vertexIndices size :"<<vertexIndices.size()<<"\n";
- // --- 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++;
- }
+ 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;
- // - 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] ---
+PView *QuadtreeLSImage::CreateLSPView(GModel* m){
- 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);
+ 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++;
- }
-
-
-
- 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++;
- }
+ }
+ std::string postTypeName = "LevelSet Value" ;
+ PView *pv = new PView (postTypeName, "NodeData", m, data, 0.0);
+ return pv;
- // --- 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();
+void QuadtreeLSImage::FillHangingNodes(){
+ _root->FillHangingNodes(_ListNodes,_ListElements,_HangingNodes);
- // --- 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";
+void Box2D::FillHangingNodes(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements,std::map< int ,std::vector<int> > &HangingNodes){
- 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::vector<std::vector<int> >::iterator it;
- return gmod;
+ 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++;
+ }
}
-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::FillMeshInfo(int & pas){
- if (!_ListNodes.empty()){
- return;
- }
+ if (!_ListNodes.empty()){
+ return;
+ }
- _root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue,pas);
+ _root->FillMeshInfo(_ListNodes,_ListElements,_ListLSValue,pas);
}
@@ -606,89 +661,88 @@ void QuadtreeLSImage::FillMeshInfo(int & 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 = _quadtree->GetRoot()->BoxSize();
-
- int xpas =pas; // Ã modifier si sizeZ n'est pas le plus petit
- int ypas = pas;
+void Box2D::FillMeshInfo(std::vector< std::vector<int> > &ListNodes, std::vector< std::vector<int> > &ListElements, std::vector<float> &ListLSValue,int &pas){
- std::cout<<"\nxyzpas :"<< xpas << " " << ypas ;
+ int size = _quadtree->GetRoot()->BoxSize();
- std::vector<std::vector<int> >::iterator it;
+ int xpas =pas; // Ã modifier si sizeZ n'est pas le plus petit
+ int ypas = pas;
- it = ListNodes.begin();
- int iti=0;
+ 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<Box*> Leafs;
- Leafs.clear();
- GetLeafsWith(Leafs, x, y);
- 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(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);
+ 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 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;
- 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);
+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;
- }
+ return;
+ }
- for (int n = 0 ; n < 4 ; n++){
- this->_children[n]->FillElementsNode(ListElements);
- }
+ for (int n = 0 ; n < 4 ; n++){
+ this->_children[n]->FillElementsNode(ListElements);
+ }
}
@@ -765,7 +819,7 @@ std::vector< std::vector <int> > *QuadtreeLSImage::GetListElements(){
/*
-void Box::GetLeafElements(std::vector< std::vector<int> > &ListElements){
+void Box2D::GetLeafElements(std::vector< std::vector<int> > &ListElements){
// if (!this->HaveChildren()) {
// std::vector<std::vector<int> >::iterator it;
@@ -802,161 +856,161 @@ void Box::GetLeafElements(std::vector< std::vector<int> > &ListElements){
void QuadtreeLSImage::SetLeafNumber(){
_LeafNumber = 0;
- _root->CountLeafNumber(_LeafNumber);
+ _root->CountLeafNumber(_LeafNumber);
}
-void Box::CountLeafNumber(int &LeafNumber){
+void Box2D::CountLeafNumber(int &LeafNumber){
- if (this->GetChild(0)==NULL) {
- LeafNumber = LeafNumber + 1;
- return;
- }
+ if (this->GetChild(0)==NULL) {
+ LeafNumber = LeafNumber + 1;
+ return;
+ }
- for (int n = 0 ; n < 4 ; n++){
- this->_children[n]->CountLeafNumber(LeafNumber);
- }
+ for (int n = 0 ; n < 4 ; n++){
+ this->_children[n]->CountLeafNumber(LeafNumber);
+ }
}
// Divide function
-void Box::Divide(){
+void Box2D::Divide(){
// box divide in 4 boxes with limits xl = (xmin + xmax)/2 ; yl = (ymin + ymax)/2
- int n;
+ int n;
- // Box 1 : < xl ; < yl
+ // Box 1 : < xl ; < yl
- n = 0;
- BoxData* data1 = new BoxData;
+ 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;
+ 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 Box(data1,this,this->_quadtree);
- this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
+ this->_children[n] = new Box2D(data1,this,this->_quadtree);
+ this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
- // Box 2 : > xl ; < yl
+ // Box 2 : > xl ; < yl
- n = 1;
- BoxData* data2 = new BoxData;
+ 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;
+ 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 Box(data2,this,this->_quadtree);
- this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
+ this->_children[n] = new Box2D(data2,this,this->_quadtree);
+ this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
- // Box 3 : < xl ; > yl
+ // Box 3 : < xl ; > yl
- n = 2;
- BoxData* data3 = new BoxData;
+ 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];
+ 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 Box(data3,this,this->_quadtree);
- this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
+ this->_children[n] = new Box2D(data3,this,this->_quadtree);
+ this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
- // Box 4 : > xl ; > yl
+ // Box 4 : > xl ; > yl
- n = 3;
- BoxData* data4 = new BoxData;
+ 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];
+ 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 Box(data4,this,this->_quadtree);
- this->_children[n]->FillLevelSetValue(_children[n]->_quadtree->GetImage(),_children[n]->GetData());
+ 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();
+ return _root->Smooth();
}
-bool Box::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;
+ }
- 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++ ){
- 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;
+ for (int n = 0 ; n < 4 ; n++){
+ if (!this->_children[n]->Smooth()) smoothed = false;
+ }
+ return smoothed;
}
// Quadtree research, box specialisation from root
-void Box::GetLeafsWith(std::vector<Box*> &Leafs, int x, int y){
+void Box2D::GetLeafsWith(std::vector<Box2D*> &Leafs, int x, int y){
- if (!this->HaveChildren()){
- Box* parent = this;
- BoxData* data;
- while(parent->GetParent()!=NULL){
- data = parent->GetParent()->GetData();
- parent=parent->GetParent();
+ if (!this->HaveChildren()){
+ Box2D* parent = this;
+ Box2DData* data;
+ while (parent->GetParent()!=NULL){
+ data = parent->GetParent()->GetData();
+ parent=parent->GetParent();
}
- Leafs.push_back(this);
- return;
- }
+ 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[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[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[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);
+ 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
index f956762c9f..d7f3eb0257 100644
--- a/contrib/arc/Classes/QuadtreeLSImage.h
+++ b/contrib/arc/Classes/QuadtreeLSImage.h
@@ -21,14 +21,14 @@
//data defining a box
-class BoxData{
+class Box2DData{
public:
int boundingbox[4]; // x1, x2, y1, y2
float LevelSetValue[2][2]; // LevelSetValue[x][y][z]
};
class QuadtreeLSImage;
-class Box;
+class Box2D;
// Quadtree class to decompose a float valued 3d image
@@ -38,7 +38,7 @@ class QuadtreeLSImage
private :
// Quadtree's root
- Box* _root;
+ Box2D* _root;
// Quadtree's image
itk::Image< float, 2 >::Pointer _image;
// Region size
@@ -49,18 +49,21 @@ class QuadtreeLSImage
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();
- Box* GetRoot(){return _root;}
+ 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;}
@@ -69,32 +72,35 @@ class QuadtreeLSImage
bool Smooth();
//std::vector< std::vector <int> > *GetListElements();
// Create GModel
- GModel* CreateGModel(bool simplex, double facx, double facy,int & sizemax,int & sizemin);
+ 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{
+class Box2D{
private :
// children pointers
- Box* _children[4];
- Box* _parent;
- BoxData* _data;
+ Box2D* _children[4];
+ Box2D* _parent;
+ Box2DData* _data;
QuadtreeLSImage* _quadtree;
int _ElementNodes[4]; // utilisé seulement par balayage
public :
- Box(BoxData* data,Box* parent,QuadtreeLSImage* octree);
- BoxData* GetData(){return _data;}
+ 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, BoxData *data);
+ 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
@@ -104,19 +110,20 @@ class Box{
// Divide by eight the box leaf
void Divide();
void PrintLevelSetValue();
- Box* GetParent(){return _parent;};
- Box* GetChild(int n){return _children[n];};
+ 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<Box*> &Leafs, int x, int y);
+ 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/XFEMInclusion.cpp b/contrib/arc/XFEMInclusion.cpp
new file mode 100644
index 0000000000..63ce7b3782
--- /dev/null
+++ b/contrib/arc/XFEMInclusion.cpp
@@ -0,0 +1,114 @@
+//
+// 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/mainImageSolver.cpp b/contrib/arc/mainImageSolver.cpp
new file mode 100644
index 0000000000..b5ce919b8d
--- /dev/null
+++ b/contrib/arc/mainImageSolver.cpp
@@ -0,0 +1,205 @@
+//
+// 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;
+}
--
GitLab