From ae3ad7d5e562d6c035b392d1625debd09254471e Mon Sep 17 00:00:00 2001
From: Boris Sedji <sedji.boris@hotmail.com>
Date: Thu, 14 Jan 2010 14:07:11 +0000
Subject: [PATCH]
---
contrib/arc/Classes/NodeEnrichedFS.cpp | 77 +-------------------
contrib/arc/Classes/NodeEnrichedFS.h | 11 +--
contrib/arc/Classes/VectorXFEMFS.h | 7 +-
contrib/arc/Classes/XFEMelasticitySolver.cpp | 31 ++++----
4 files changed, 23 insertions(+), 103 deletions(-)
diff --git a/contrib/arc/Classes/NodeEnrichedFS.cpp b/contrib/arc/Classes/NodeEnrichedFS.cpp
index 44bd73527f..166082489a 100644
--- a/contrib/arc/Classes/NodeEnrichedFS.cpp
+++ b/contrib/arc/Classes/NodeEnrichedFS.cpp
@@ -11,31 +11,6 @@
#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)
{
@@ -63,7 +38,7 @@ template <class T> void NodeEnrichedFS<T>::f(MElement *ele, double u, double v,
}
int curpos=vals.size();
- vals.resize(curpos+nbdofs);
+ vals.reserve(curpos+nbdofs);
// first normal dofs
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,
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)
@@ -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)));
- }
- }
-};
diff --git a/contrib/arc/Classes/NodeEnrichedFS.h b/contrib/arc/Classes/NodeEnrichedFS.h
index 5767e94ae7..8eab580750 100644
--- a/contrib/arc/Classes/NodeEnrichedFS.h
+++ b/contrib/arc/Classes/NodeEnrichedFS.h
@@ -20,6 +20,7 @@
#include "Numeric.h"
#include "MElement.h"
#include "MVertex.h"
+#include "MPoint.h"
#include "dofManager.h"
#include "functionSpace.h"
#include "simpleFunction.h"
@@ -34,9 +35,7 @@ public:
typedef typename TensorialTraits<T>::ValType ValType;
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:
@@ -57,18 +56,12 @@ public:
virtual ~NodeEnrichedFS() {}
// 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
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
- virtual void getKeys(MVertex *ver, std::vector<Dof> &keys);
};
#endif
diff --git a/contrib/arc/Classes/VectorXFEMFS.h b/contrib/arc/Classes/VectorXFEMFS.h
index 23c577f38d..8a156bb5c8 100644
--- a/contrib/arc/Classes/VectorXFEMFS.h
+++ b/contrib/arc/Classes/VectorXFEMFS.h
@@ -28,9 +28,6 @@ class ScalarLagrangeToXfemFS : public ScalarLagrangeFunctionSpace{
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:
@@ -74,9 +71,7 @@ class ScalarXFEMToVectorFS : public ScalarToAnyFunctionSpace<SVector3>
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::set<int > & TagEnrichedVertex , simpleFunction<double> * funcEnrichment) :
ScalarToAnyFunctionSpace<SVector3>::ScalarToAnyFunctionSpace(ScalarLagrangeToXfemFS(id,TagEnrichedVertex,funcEnrichment),
diff --git a/contrib/arc/Classes/XFEMelasticitySolver.cpp b/contrib/arc/Classes/XFEMelasticitySolver.cpp
index 63e94e0b08..c405812926 100644
--- a/contrib/arc/Classes/XFEMelasticitySolver.cpp
+++ b/contrib/arc/Classes/XFEMelasticitySolver.cpp
@@ -26,11 +26,12 @@
#if defined(HAVE_POST)
#include "PView.h"
#include "PViewData.h"
+#include "MPoint.h"
#endif
#include "DILevelset.h"
-#include "VectorXFEMFS.cpp"
+//#include "VectorXFEMFS.cpp"
#include "XFEMelasticitySolver.h"
#include "FuncHeaviside.h"
@@ -84,7 +85,6 @@ void XFEMelasticitySolver::solve(){
_EnrichComp.push_back(1);
//_EnrichComp.push_back(2);
-
// level set definition (in .dat ??)
double a(0.), b(1.), c(0.), d(-4.7);
int n(1); // pour level set
@@ -105,13 +105,11 @@ void XFEMelasticitySolver::solve(){
// give priority to fixations : when a dof is fixed, it cannot be
// numbered afterwards
+ std::cout << "Dirichlet BC"<< std::endl;
for (unsigned int i = 0; i < allDirichlet.size(); i++)
{
FilterDofComponent filter(allDirichlet[i]._comp);
- 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);
}
// we number the dofs : when a dof is numbered, it cannot be numbered
@@ -125,22 +123,20 @@ void XFEMelasticitySolver::solve(){
// 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);
- 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);
}
- // bulk material law
+// 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);
}
@@ -158,12 +154,15 @@ void XFEMelasticitySolver::solve(){
Assemble(Elastic_Energy_Term,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,energ);
}
printf("elastic energy=%f\n",energ);
+ for (int i=0;i<pAssembler->sizeOfR();i++) std::cout<< lsys->getFromSolution(i) << "\n";
+
}
PView* XFEMelasticitySolver::buildDisplacementView (const std::string &postFileName)
{
+
std::set<MVertex*> v;
for (unsigned int i = 0; i < elasticFields.size(); ++i)
{
@@ -179,10 +178,14 @@ PView* XFEMelasticitySolver::buildDisplacementView (const std::string &postFileN
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;
+ 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;
+
}
--
GitLab