Skip to content
Snippets Groups Projects
Commit cb743a29 authored by Boris Sedji's avatar Boris Sedji
Browse files

No commit message

No commit message
parent cf2974e2
No related branches found
No related tags found
No related merge requests found
//
// Description : Heaviside function based on level set discontinuity description
//
//
// Author: <Boris Sedji>, 01/2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#ifndef _FUNCHEAVISIDE_H_
#define _FUNCHEAVISIDE_H_
#include "simpleFunction.h"
#include "DILevelset.h"
class FuncHeaviside : public simpleFunction<double> {
private :
gLevelset *_ls;
public :
FuncHeaviside(gLevelset *ls){
_ls = ls;
}
virtual double operator () (double x,double y,double z) const {
if (_ls->isInsideDomain (x,y,z)){
return 1 ;
}else{
return -1;
}
}
virtual void gradient (double x, double y, double z,
double & dfdx, double & dfdy, double & dfdz) const
{ dfdx = dfdy = dfdz = 0.0; }
};
#endif
//
// Description : XFEM elasticity solver implementation
//
//
// Author: <Eric Bechet>::<Boris Sedji>, 01/2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "VectorXFEMFS.h"
#include "SPoint3.h"
void ScalarLagrangeToXfemFS::f(MVertex *ver, std::vector<ValType> &vals)
{
if (!(*_TagEnrichedVertex)[ver->getNum()]) vals.push_back(1.0);
else
{
double func = (*_funcEnrichment)(ver->x(),ver->y(),ver->z());
vals.push_back(1.0);
vals.push_back(func);
}
}
void ScalarLagrangeToXfemFS::f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
// Enrichment function calcul
SPoint3 p;
ele->pnt(u, v, w, p);
double func;
func = (*_funcEnrichment)(p.x(), p.y(), p.z());
// We need parent parameters
if (ele->getParent()) ele = ele->getParent();
int ndofs,normaldofs;
ndofs = ele->getNumVertices();
normaldofs = ndofs;
for (int i=0 ;i<ele->getNumVertices();i++)
{
if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()]) ndofs = ndofs + 1; // enriched dof
}
int curpos=vals.size();
vals.resize(curpos+ndofs);
// normal shape functions
ele->getShapeFunctions(u, v, w, &(vals[curpos]));
int k=0;
for (int i=0 ;i<ele->getNumVertices();i++)
{
if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()])
{
vals[curpos+normaldofs+k] = vals[curpos+i]*func; // enriched dof
k++;
}
}
};
void ScalarLagrangeToXfemFS::gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
{
// Enrichment function calcul
SPoint3 p;
ele->pnt(u, v, w, p);
double func;
func = (*_funcEnrichment)(p.x(), p.y(), p.z());
// We need parent parameters
if (ele->getParent()) ele = ele->getParent();
int ndofs,normaldofs;
ndofs= ele->getNumVertices();
normaldofs = ndofs; // normal dofs
for (int i=0 ;i< (ele->getNumVertices());i++)
{
if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()]) ndofs = ndofs + 1; // enriched dof
}
int curpos = grads.size();
grads.reserve(curpos+ndofs);
double gradsuvw[256][3];
ele->getGradShapeFunctions(u, v, w, gradsuvw);
int k = 0;
for (int i=0 ;i<(ele->getNumVertices());i++)
{
if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()]) // enriched dof
{
gradsuvw[curpos+normaldofs+k][0] = gradsuvw[curpos+i][0]*func;
gradsuvw[curpos+normaldofs+k][1] = gradsuvw[curpos+i][1]*func;
gradsuvw[curpos+normaldofs+k][2] = gradsuvw[curpos+i][2]*func;
k++;
}
}
double jac[3][3];
double invjac[3][3];
const double detJ = ele->getJacobian(u, v, w, jac);
// redondant : on fait cet appel a l'exterieur
inv3x3(jac, invjac);
for (int i=0;i<ndofs;++i)
grads.push_back(GradType(
invjac[0][0] * gradsuvw[i][0] + invjac[0][1] * gradsuvw[i][1] + invjac[0][2] * gradsuvw[i][2],
invjac[1][0] * gradsuvw[i][0] + invjac[1][1] * gradsuvw[i][1] + invjac[1][2] * gradsuvw[i][2],
invjac[2][0] * gradsuvw[i][0] + invjac[2][1] * gradsuvw[i][1] + invjac[2][2] * gradsuvw[i][2]
));
};
int ScalarLagrangeToXfemFS::getNumKeys(MElement *ele)
{
// We need parent parameters
if (ele->getParent()) ele = ele->getParent();
int ndofs;
ndofs = ele->getNumVertices();
{
for (int i=0 ;i<(ele->getNumVertices());i++)
{
if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()]) ndofs = ndofs + 1;
}
}
return ndofs;
}
int ScalarLagrangeToXfemFS::getNumKeys(MVertex *ver)
{
if (!(*_TagEnrichedVertex)[ver->getNum()]) return 1 ;
else return 2; // if enriched vertex, there is two dof at this vertex (one dimension)
}
void ScalarLagrangeToXfemFS::getKeys(MElement *ele, std::vector<Dof> &keys)
{
if (ele->getParent()) ele = ele->getParent();
// vector of additionals keys
std::vector<Dof> SupKeys;
// vector of all vertex dof keys
std::vector<Dof> VertexKeys;
int ndofs,normaldofs;
ndofs = ele->getNumVertices();
for (int i=0 ;i<(ele->getNumVertices());i++)
{
if ((*_TagEnrichedVertex)[ele->getVertex(i)->getNum()]) ndofs = ndofs + 1;
}
keys.reserve(keys.size()+ndofs);
// the order of dof in one element of three vertex is u1x u2x u3x a1x a2x a3x u1y, etc...
for (int i=0;i< (ele->getNumVertices()) ;++i)
{
getKeys(ele->getVertex(i), VertexKeys);
if (VertexKeys.size()>1) SupKeys.push_back(VertexKeys[1]);
keys.push_back(VertexKeys[0]);
VertexKeys.clear();
}
for (int i=0;i< SupKeys.size() ;++i) keys.push_back(SupKeys[i]);
}
void ScalarLagrangeToXfemFS::getKeys(MVertex *ver, std::vector<Dof> &keys)
{
if ((*_TagEnrichedVertex)[ver->getNum()])
{
keys.push_back(Dof(ver->getNum(), _iField));
keys.push_back(Dof(ver->getNum(), _iField+1)); // we tag the additional dof with number field+1
}
else keys.push_back(Dof(ver->getNum(), _iField));
};
//
// Description : From Scalar Lagrange Function Space To XFEM on Tagged Vertex
//
//
// Author: <Eric Bechet>::<Boris Sedji>, 01/2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "SVector3.h"
#include "STensor3.h"
#include <vector>
#include <iterator>
#include <iostream>
#include "Numeric.h"
#include "MElement.h"
#include "dofManager.h"
#include "functionSpace.h"
#include "simpleFunction.h"
class ScalarLagrangeToXfemFS : public ScalarLagrangeFunctionSpace{
public:
typedef TensorialTraits<double>::ValType ValType;
typedef TensorialTraits<double>::GradType GradType;
typedef TensorialTraits<double>::HessType HessType;
typedef TensorialTraits<double>::DivType DivType;
typedef TensorialTraits<double>::CurlType CurlType;
protected:
std::map<int , bool> *_TagEnrichedVertex;
simpleFunction<double> *_funcEnrichment;
public:
//
ScalarLagrangeToXfemFS(int i, std::map<int , bool> &TagEnrichedVertex, simpleFunction<double> *funcEnrichment) : ScalarLagrangeFunctionSpace(i)
{
_TagEnrichedVertex = & TagEnrichedVertex;
_funcEnrichment = funcEnrichment;
}
virtual int getId(void) const {return ScalarLagrangeFunctionSpace::_iField;};
// Shape function value in element at uvw
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals);
// Shape function value at vertex
virtual void f(MVertex *ver, std::vector<ValType> &vals) ;
// Grad Shape function value in element at uvw
virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) ;
// Shape function Number for element (in one dimension (scalar))
virtual int getNumKeys(MElement *ele);
// Dof keys for the element
virtual void getKeys(MElement *ele, std::vector<Dof> &keys) ;
// Shape function Number associate with the vertex
virtual int getNumKeys(MVertex *ver);
// Get Dof keys for the vertex (in one dimension (scalar))
virtual void getKeys(MVertex *ver, std::vector<Dof> &keys);
};
class ScalarXFEMToVectorFS : public ScalarToAnyFunctionSpace<SVector3>
{
protected:
static const SVector3 BasisVectors[3];
public:
enum Along { VECTOR_X=0, VECTOR_Y=1, VECTOR_Z=2 };
typedef TensorialTraits<SVector3>::ValType ValType;
typedef TensorialTraits<SVector3>::GradType GradType;
typedef TensorialTraits<SVector3>::HessType HessType;
typedef TensorialTraits<SVector3>::DivType DivType;
typedef TensorialTraits<SVector3>::CurlType CurlType;
ScalarXFEMToVectorFS(int id , std::map<int , bool> & TagEnrichedVertex , simpleFunction<double> * funcEnrichment) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment),
SVector3(1.,0.,0.),VECTOR_X, SVector3(0.,1.,0.),VECTOR_Y,SVector3(0.,0.,1.),VECTOR_Z) {}
ScalarXFEMToVectorFS(int id,Along comp1, std::map<int , bool> &TagEnrichedVertex , simpleFunction<double> * funcEnrichment) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment),
BasisVectors[comp1],comp1) {}
ScalarXFEMToVectorFS(int id,Along comp1,Along comp2, std::map<int , bool> &TagEnrichedVertex,simpleFunction<double> *funcEnrichment) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment),
BasisVectors[comp1],comp1, BasisVectors[comp2],comp2) {}
ScalarXFEMToVectorFS(int id,Along comp1,Along comp2, Along comp3, std::map<int , bool> &TagEnrichedVertex,simpleFunction<double> *funcEnrichment) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment),
BasisVectors[comp1],comp1, BasisVectors[comp2],comp2, BasisVectors[comp3],comp3) {}
};
const SVector3 ScalarXFEMToVectorFS::BasisVectors[3]={SVector3(1,0,0),SVector3(0,1,0),SVector3(0,0,1)};
//
// Description : XFEM elasticity solver, element space function enriched on tagged vertex
//
//
// Author: <Eric Bechet>::<Boris Sedji>, 01/2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include <string.h>
#include "GmshConfig.h"
#include "elasticitySolver.h"
#include "linearSystemCSR.h"
#include "linearSystemPETSc.h"
#include "linearSystemGMM.h"
#include "Numeric.h"
#include "GModel.h"
#include "functionSpace.h"
#include "terms.h"
#include "solverAlgorithms.h"
#include "quadratureRules.h"
#include "solverField.h"
#if defined(HAVE_POST)
#include "PView.h"
#include "PViewData.h"
#endif
#include "DILevelset.h"
#include "VectorXFEMFS.cpp"
#include "XFEMelasticitySolver.h"
#include "FuncHeaviside.h"
XFEMelasticitySolver::~XFEMelasticitySolver()
{
delete _funcEnrichment;
}
void XFEMelasticitySolver::setMesh(const std::string &meshFileName)
{
pModel = new GModel();
pModel->readMSH(meshFileName.c_str());
_dim = pModel->getNumRegions() ? 3 : 2;
}
void XFEMelasticitySolver::solve(){
#if defined(HAVE_TAUCS)
linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
#elif defined(HAVE_PETSC)
linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
#else
linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
lsys->setNoisy(2);
#endif
if (pAssembler) delete pAssembler;
pAssembler = new dofManager<double>(lsys);
// determine all enriched nodes and save in map member _TagEnrichedVertex
for (int i = 0; i < elasticFields.size(); ++i){
groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
for ( ; it != elasticFields[i].g->end(); ++it)
{
MElement *e = *it;
if (e->getParent()) // if element got parents
{
for (int k = 0; k < e->getParent()->getNumVertices(); ++k)
{ // for all vertices in the element parent
_TagEnrichedVertex[e->getParent()->getVertex(k)->getNum()] = true;
}
}
else
{
for (int k = 0; k < e->getNumVertices(); ++k)
{
_TagEnrichedVertex[e->getVertex(k)->getNum()] = false;
}
}
}
}
// level set definition (in .dat ??)
double a(0.), b(1.), c(0.), d(-4.7);
int n(1); // pour level set
gLevelset *ls = new gLevelsetPlane(a, b, c, d, n);
_funcEnrichment = new FuncHeaviside(ls);
// space function definition
if (LagSpace) delete LagSpace;
if (_dim==3) LagSpace=new ScalarXFEMToVectorFS(_tag,_TagEnrichedVertex,_funcEnrichment);
if (_dim==2) LagSpace=new ScalarXFEMToVectorFS(_tag,ScalarXFEMToVectorFS::VECTOR_X,ScalarXFEMToVectorFS::VECTOR_Y,_TagEnrichedVertex,_funcEnrichment);
// we first do all fixations. the behavior of the dofManager is to
// give priority to fixations : when a dof is fixed, it cannot be
// numbered afterwards
for (unsigned int i = 0; i < allDirichlet.size(); i++)
{
FilterDofComponent filter(allDirichlet[i]._comp);
if (allDirichlet[i].onWhat==BoundaryCondition::ON_VERTEX)
FixNodalDofs(*LagSpace,allDirichlet[i].g->vbegin(),allDirichlet[i].g->vend(),*pAssembler,allDirichlet[i]._f,filter);
else
FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,allDirichlet[i]._f,filter);
}
// we number the dofs : when a dof is numbered, it cannot be numbered
// again with another number.
for (unsigned int i = 0; i < elasticFields.size(); ++i)
{
NumberDofs(*LagSpace, elasticFields[i].g->begin(), elasticFields[i].g->end(),*pAssembler);
}
// Now we start the assembly process
// First build the force vector
GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
for (unsigned int i = 0; i < allNeumann.size(); i++)
{
LoadTerm<SVector3> Lterm(*LagSpace,allNeumann[i]._f);
if (allNeumann[i].onWhat==BoundaryCondition::ON_VERTEX)
Assemble(Lterm,*LagSpace,allNeumann[i].g->vbegin(),allNeumann[i].g->vend(),*pAssembler);
else
Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler);
}
// bulk material law
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); i++)
{
IsotropicElasticTerm Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu);
Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
}
printf("-- done assembling!\n");
std::cout<<"Dof Number : " << pAssembler->sizeOfR() <<"\n";
lsys->systemSolve();
printf("-- done solving!\n");
double energ=0;
for (unsigned int i = 0; i < elasticFields.size(); i++)
{
SolverField<SVector3> Field(pAssembler, LagSpace);
IsotropicElasticTerm Eterm(Field,elasticFields[i]._E,elasticFields[i]._nu);
BilinearTermToScalarTerm<SVector3,SVector3> Elastic_Energy_Term(Eterm);
Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
}
printf("elastic energy=%f\n",energ);
}
PView* XFEMelasticitySolver::buildDisplacementView (const std::string &postFileName)
{
std::set<MVertex*> v;
for (unsigned int i = 0; i < elasticFields.size(); ++i)
{
for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
{
MElement *e=*it;
if (e->getParent()) e=e->getParent();
for (int j = 0; j < e->getNumVertices(); ++j) v.insert(e->getVertex(j));
}
}
std::map<int, std::vector<double> > data;
SolverField<SVector3> Field(pAssembler, LagSpace);
for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
{
SVector3 val;
Field.f(*it,val);
std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=0;
data[(*it)->getNum()]=vec;
}
PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
return pv;
}
//
// Description : XFEM elasticity solver, element space function enriched on tagged vertex
//
//
// Author: <Eric Bechet>::<Boris Sedji>, 01/2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#ifndef _XFEMELASTICITY_SOLVER_H_
#define _XFEMELASTICITY_SOLVER_H_
#include "elasticitySolver.h"
#include "simpleFunction.h"
class XFEMelasticitySolver : public elasticitySolver
{
protected :
// map containing the tag of vertex and enriched status
std::map<int , bool> _TagEnrichedVertex;
// simple multiplying function enrichment to enrich the space function
simpleFunction <double> *_funcEnrichment;
public :
XFEMelasticitySolver(int tag) : elasticitySolver(tag) {}
~XFEMelasticitySolver();
// create a GModel and determine de dimension of mesh in meshFileName
virtual void setMesh(const std::string &meshFileName);
// system solve, read the .dat file, fill tagEnrichedVertex, fill funcEnrichment, solve
virtual void solve();
virtual PView *buildDisplacementView(const std::string &postFileName);
};
#endif
//
// Description : Import Float Level Set Image in Gmsh QuadTree Mesh
//
//
// Author: <Boris Sedji>, 12/2009
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include <stdio.h>
#include "GModel.h"
#include "Classes/QuadtreeLSImage.cpp"
int main( int argc, char *argv[] )
{
if( argc < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " inputLevelSet2dImage MeshSizeMax MeshSizeMin"<< std::endl;
return 1;
}
typedef float PixelTypeFloat;
const unsigned int Dimension = 2;
typedef itk::Image< PixelTypeFloat, Dimension > ImageTypeFloat;
typedef itk::ImageFileReader< ImageTypeFloat > ImageReaderTypeFloat;
ImageReaderTypeFloat::Pointer reader = ImageReaderTypeFloat::New();
reader->SetFileName(argv[1]);
ImageTypeFloat::Pointer image = reader->GetOutput();
image->Update();
ImageTypeFloat::RegionType region;
region = image->GetLargestPossibleRegion ();
std::cout<<"\nImage dimensions : " << region.GetSize(0) << " x " << region.GetSize(1);
QuadtreeLSImage quadtree(image);
int sizemax = atoi(argv[2]);
int sizemin = atoi(argv[3]);
// Mesh with conditions on mesh size
quadtree.Mesh(sizemax,sizemin);
// Smoothing mesh to avoid too much generation between adjacent elements
bool statut=false;
int k = 0;
std::cout<<"\nLeaf Number : "<< (quadtree.GetLeafNumber())<<"\n";
while((!statut) & (k < 20)){
statut = quadtree.Smooth();
std::cout<<"\nk : "<< k;
std::cout<<"\nsmoothed : " << (int)statut;
std::cout<<"\nLeaf Number : "<< (quadtree.GetLeafNumber())<<"\n";
k++;
}
// Create GModel with the octree mesh
GModel *m = quadtree.CreateGModel();
// Write a .msh file with the mesh
std::string ModelName = "QuadtreeMesh.msh" ;
m->writeMSH(ModelName,2.1,false,false);
// Write a .msh file with the level set values as postview data
PView *pv = quadtree.CreateLSPView(m);
bool useadapt = true;
pv->getData(useadapt)->writeMSH("LSPView.msh", false);
std::cout<<"\n";
return 0;
}
//
// Description : XFEM elasticity solver main
//
//
// Author: <Eric Bechet>::<Boris Sedji>, 01/2010
//
// Copyright: See COPYING file that comes with this distribution
//
//
#include "Gmsh.h"
#include "elasticitySolver.h"
#include "PView.h"
#include "PViewData.h"
#include "highlevel.h"
#include "groupOfElements.h"
#include <iterator>
#include "Classes/XFEMelasticitySolver.cpp"
int main (int argc, char* argv[])
{
if (argc < 2)
{
printf("usage : elasticity input_file_name yeahhh \n");
return -1;
}
GmshInitialize(argc, argv);
// globals are still present in Gmsh
// instanciate a solver
XFEMelasticitySolver mySolver (1000);
// read some input file
mySolver.readInputFile(argv[1]);
// solve the problem
mySolver.solve();
PView *pv = mySolver.buildDisplacementView("displacement");
pv->getData()->writeMSH("disp.msh", false);
delete pv;
pv = mySolver.buildElasticEnergyView("elastic energy");
pv->getData()->writeMSH("energ.msh", false);
delete pv;
// stop gmsh
GmshFinalize();
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment