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

No commit message

No commit message
parent 9036fc14
No related branches found
No related tags found
No related merge requests found
...@@ -11,31 +11,6 @@ ...@@ -11,31 +11,6 @@
#include "NodeEnrichedFS.h" #include "NodeEnrichedFS.h"
template <class T> void NodeEnrichedFS<T>::f(MVertex *ver, std::vector<ValType> &vals)
{
std::vector<ValType> valsd;
// get the support value at the vertex
_SupportFS->f(ver,valsd);
int normaldofs=valsd.size();
int curpos = vals.size();
std::set<int>::iterator it;
it = _TagEnrichedVertex->find(ver->getNum());
if (it!=_TagEnrichedVertex->end()) // if node enriched
{
double func = (*_funcEnrichment)(ver->x(),ver->y(),ver->z());
int nbdofs = normaldofs + _EnrichComp->size();
vals.reserve(curpos+nbdofs);
for (int i=0;i<normaldofs;i++) vals.push_back(valsd[i]);
for (int i=0;i<_EnrichComp->size();i++) vals.push_back(valsd[_EnrichComp->at(i)]*func);
}
else
{
int nbdofs = normaldofs;
vals.reserve(curpos+nbdofs);
for (int i=0;i<normaldofs;i++) vals.push_back(valsd[i]);
}
}
template <class T> void NodeEnrichedFS<T>::f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals) template <class T> void NodeEnrichedFS<T>::f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{ {
...@@ -63,7 +38,7 @@ template <class T> void NodeEnrichedFS<T>::f(MElement *ele, double u, double v, ...@@ -63,7 +38,7 @@ template <class T> void NodeEnrichedFS<T>::f(MElement *ele, double u, double v,
} }
int curpos=vals.size(); int curpos=vals.size();
vals.resize(curpos+nbdofs); vals.reserve(curpos+nbdofs);
// first normal dofs // first normal dofs
for (int i=0;i<normaldofs;i++) vals.push_back(valsd[i]); for (int i=0;i<normaldofs;i++) vals.push_back(valsd[i]);
...@@ -82,6 +57,8 @@ template <class T> void NodeEnrichedFS<T>::f(MElement *ele, double u, double v, ...@@ -82,6 +57,8 @@ template <class T> void NodeEnrichedFS<T>::f(MElement *ele, double u, double v,
for (int j=0;j<_EnrichComp->size();j++) vals.push_back(valsd[EnrichedVertex[i]+(_EnrichComp->at(j))*elep->getNumVertices()]*func); for (int j=0;j<_EnrichComp->size();j++) vals.push_back(valsd[EnrichedVertex[i]+(_EnrichComp->at(j))*elep->getNumVertices()]*func);
} }
} }
} }
template <class T> void NodeEnrichedFS<T>::gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) template <class T> void NodeEnrichedFS<T>::gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
...@@ -204,51 +181,3 @@ template <class T> void NodeEnrichedFS<T>::getKeys(MElement *ele, std::vector<Do ...@@ -204,51 +181,3 @@ template <class T> void NodeEnrichedFS<T>::getKeys(MElement *ele, std::vector<Do
} }
} }
} }
template <class T> int NodeEnrichedFS<T>::getNumKeys(MVertex *ver)
{
std::set<int>::iterator it;
it = _TagEnrichedVertex->find(ver->getNum());
if (it!=_TagEnrichedVertex->end()) // if node enriched
{
return _SupportFS->getNumKeys(ver)+_EnrichComp->size();
}
else
{
return _SupportFS->getNumKeys(ver);
}
}
template <class T> void NodeEnrichedFS<T>::getKeys(MVertex *ver, std::vector<Dof> &keys)
{
int normalkeys= _SupportFS->getNumKeys(ver);
std::vector<Dof> bufk;
bufk.reserve(normalkeys);
_SupportFS->getKeys(ver,bufk);
int normaldofs=bufk.size();
int nbdofs = normaldofs;
int curpos=keys.size();
std::set<int>::iterator it;
it = _TagEnrichedVertex->find(ver->getNum());
if (it!=_TagEnrichedVertex->end())
nbdofs = nbdofs+_EnrichComp->size(); // enriched dof
keys.reserve(curpos+nbdofs);
// normal dofs
for (int i=0;i<bufk.size();i++) keys.push_back(bufk[i]);
// enriched dofs
if (nbdofs>normaldofs)
{
for (int j=0;j<_EnrichComp->size();j++)
{
int i1,i2;
Dof::getTwoIntsFromType(bufk[_EnrichComp->at(j)].getType(), i1,i2);
keys.push_back(Dof(bufk[_EnrichComp->at(j)].getEntity(),Dof::createTypeWithTwoInts(_EnrichComp->at(j),i1+1)));
}
}
};
...@@ -20,6 +20,7 @@ ...@@ -20,6 +20,7 @@
#include "Numeric.h" #include "Numeric.h"
#include "MElement.h" #include "MElement.h"
#include "MVertex.h" #include "MVertex.h"
#include "MPoint.h"
#include "dofManager.h" #include "dofManager.h"
#include "functionSpace.h" #include "functionSpace.h"
#include "simpleFunction.h" #include "simpleFunction.h"
...@@ -34,9 +35,7 @@ public: ...@@ -34,9 +35,7 @@ public:
typedef typename TensorialTraits<T>::ValType ValType; typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType; typedef typename TensorialTraits<T>::GradType GradType;
typedef typename TensorialTraits<T>::HessType HessType;
typedef typename TensorialTraits<T>::DivType DivType;
typedef typename TensorialTraits<T>::CurlType CurlType;
protected: protected:
...@@ -57,18 +56,12 @@ public: ...@@ -57,18 +56,12 @@ public:
virtual ~NodeEnrichedFS() {} virtual ~NodeEnrichedFS() {}
// Shape function value in element at uvw // Shape function value in element at uvw
virtual void f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals); 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 // Grad Shape function value in element at uvw
virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) ; virtual void gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads) ;
// Shape function number for element // Shape function number for element
virtual int getNumKeys(MElement *ele); virtual int getNumKeys(MElement *ele);
// Dof keys for the element // Dof keys for the element
virtual void getKeys(MElement *ele, std::vector<Dof> &keys) ; 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
virtual void getKeys(MVertex *ver, std::vector<Dof> &keys);
}; };
#endif #endif
...@@ -28,9 +28,6 @@ class ScalarLagrangeToXfemFS : public ScalarLagrangeFunctionSpace{ ...@@ -28,9 +28,6 @@ class ScalarLagrangeToXfemFS : public ScalarLagrangeFunctionSpace{
typedef TensorialTraits<double>::ValType ValType; typedef TensorialTraits<double>::ValType ValType;
typedef TensorialTraits<double>::GradType GradType; typedef TensorialTraits<double>::GradType GradType;
typedef TensorialTraits<double>::HessType HessType;
typedef TensorialTraits<double>::DivType DivType;
typedef TensorialTraits<double>::CurlType CurlType;
protected: protected:
...@@ -74,9 +71,7 @@ class ScalarXFEMToVectorFS : public ScalarToAnyFunctionSpace<SVector3> ...@@ -74,9 +71,7 @@ class ScalarXFEMToVectorFS : public ScalarToAnyFunctionSpace<SVector3>
typedef TensorialTraits<SVector3>::ValType ValType; typedef TensorialTraits<SVector3>::ValType ValType;
typedef TensorialTraits<SVector3>::GradType GradType; typedef TensorialTraits<SVector3>::GradType GradType;
typedef TensorialTraits<SVector3>::HessType HessType;
typedef TensorialTraits<SVector3>::DivType DivType;
typedef TensorialTraits<SVector3>::CurlType CurlType;
ScalarXFEMToVectorFS(int id , std::set<int > & TagEnrichedVertex , simpleFunction<double> * funcEnrichment) : ScalarXFEMToVectorFS(int id , std::set<int > & TagEnrichedVertex , simpleFunction<double> * funcEnrichment) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment), ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment),
......
...@@ -26,11 +26,12 @@ ...@@ -26,11 +26,12 @@
#if defined(HAVE_POST) #if defined(HAVE_POST)
#include "PView.h" #include "PView.h"
#include "PViewData.h" #include "PViewData.h"
#include "MPoint.h"
#endif #endif
#include "DILevelset.h" #include "DILevelset.h"
#include "VectorXFEMFS.cpp" //#include "VectorXFEMFS.cpp"
#include "XFEMelasticitySolver.h" #include "XFEMelasticitySolver.h"
#include "FuncHeaviside.h" #include "FuncHeaviside.h"
...@@ -84,7 +85,6 @@ void XFEMelasticitySolver::solve(){ ...@@ -84,7 +85,6 @@ void XFEMelasticitySolver::solve(){
_EnrichComp.push_back(1); _EnrichComp.push_back(1);
//_EnrichComp.push_back(2); //_EnrichComp.push_back(2);
// level set definition (in .dat ??) // level set definition (in .dat ??)
double a(0.), b(1.), c(0.), d(-4.7); double a(0.), b(1.), c(0.), d(-4.7);
int n(1); // pour level set int n(1); // pour level set
...@@ -105,12 +105,10 @@ void XFEMelasticitySolver::solve(){ ...@@ -105,12 +105,10 @@ void XFEMelasticitySolver::solve(){
// give priority to fixations : when a dof is fixed, it cannot be // give priority to fixations : when a dof is fixed, it cannot be
// numbered afterwards // numbered afterwards
std::cout << "Dirichlet BC"<< std::endl;
for (unsigned int i = 0; i < allDirichlet.size(); i++) for (unsigned int i = 0; i < allDirichlet.size(); i++)
{ {
FilterDofComponent filter(allDirichlet[i]._comp); 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); FixNodalDofs(*LagSpace,allDirichlet[i].g->begin(),allDirichlet[i].g->end(),*pAssembler,allDirichlet[i]._f,filter);
} }
...@@ -125,13 +123,10 @@ void XFEMelasticitySolver::solve(){ ...@@ -125,13 +123,10 @@ void XFEMelasticitySolver::solve(){
// First build the force vector // First build the force vector
GaussQuadrature Integ_Boundary(GaussQuadrature::Val); GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
std::cout << "Neumann BC"<< std::endl;
for (unsigned int i = 0; i < allNeumann.size(); i++) for (unsigned int i = 0; i < allNeumann.size(); i++)
{ {
LoadTerm<SVector3> Lterm(*LagSpace,allNeumann[i]._f); 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); Assemble(Lterm,*LagSpace,allNeumann[i].g->begin(),allNeumann[i].g->end(),Integ_Boundary,*pAssembler);
} }
...@@ -141,6 +136,7 @@ void XFEMelasticitySolver::solve(){ ...@@ -141,6 +136,7 @@ void XFEMelasticitySolver::solve(){
for (unsigned int i = 0; i < elasticFields.size(); i++) for (unsigned int i = 0; i < elasticFields.size(); i++)
{ {
IsotropicElasticTerm Eterm(*LagSpace,elasticFields[i]._E,elasticFields[i]._nu); 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); Assemble(Eterm,*LagSpace,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
} }
...@@ -158,12 +154,15 @@ void XFEMelasticitySolver::solve(){ ...@@ -158,12 +154,15 @@ void XFEMelasticitySolver::solve(){
Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ); Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
} }
printf("elastic energy=%f\n",energ); printf("elastic energy=%f\n",energ);
for (int i=0;i<pAssembler->sizeOfR();i++) std::cout<< lsys->getFromSolution(i) << "\n";
} }
PView* XFEMelasticitySolver::buildDisplacementView (const std::string &postFileName) PView* XFEMelasticitySolver::buildDisplacementView (const std::string &postFileName)
{ {
std::set<MVertex*> v; std::set<MVertex*> v;
for (unsigned int i = 0; i < elasticFields.size(); ++i) for (unsigned int i = 0; i < elasticFields.size(); ++i)
{ {
...@@ -179,10 +178,14 @@ PView* XFEMelasticitySolver::buildDisplacementView (const std::string &postFileN ...@@ -179,10 +178,14 @@ PView* XFEMelasticitySolver::buildDisplacementView (const std::string &postFileN
for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it) for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
{ {
SVector3 val; SVector3 val;
Field.f(*it,val); MPoint p(*it);
std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=0; 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; data[(*it)->getNum()]=vec;
} }
PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0); PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
return pv; return pv;
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment