From e546b03d89b7ec06f28fd2dd28a700d8ff81290c Mon Sep 17 00:00:00 2001
From: Eric Bechet <eric.bechet@ulg.ac.be>
Date: Sat, 12 Dec 2009 08:21:24 +0000
Subject: [PATCH] -added assembly for linear terms in dofmanager -done BC with
linear terms in elasiticiy
---
Solver/dofManager.h | 81 +++++++++++-------
Solver/elasticitySolver.cpp | 160 ++++++++++++++++++++++++++----------
Solver/functionSpace.h | 2 +-
Solver/terms.h | 38 ++++++++-
4 files changed, 209 insertions(+), 72 deletions(-)
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 5ffcbf3b69..91bd3cf76b 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -179,48 +179,73 @@ class dofManager{
}
for (unsigned int i = 0; i < R.size(); i++){
if (NR[i] != -1){
- for (unsigned int j = 0; j < C.size(); j++){
- if (NC[j] != -1){
- _current->addToMatrix(NR[i], NC[j], m(i, j));
- }
- else{
- typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C[j]);
- if (itFixed != fixed.end()){
- _current->addToRightHandSide(NR[i], -m(i,j) * itFixed->second);
- }
- }
- }
+ for (unsigned int j = 0; j < C.size(); j++){
+ if (NC[j] != -1){
+ _current->addToMatrix(NR[i], NC[j], m(i, j));
+ }
+ else{
+ typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(C[j]);
+ if (itFixed != fixed.end()){
+ _current->addToRightHandSide(NR[i], -m(i,j) * itFixed->second);
+ }
+ }
+ }
}
}
}
- inline void assemble(std::vector<Dof> &R, fullMatrix<dataMat> &m)
+
+ inline void assemble(std::vector<Dof> &R, fullVector<dataMat> &m) // for linear forms
{
if (!_current->isAllocated()) _current->allocate(unknown.size());
-
std::vector<int> NR(R.size());
-
- for (unsigned int i = 0; i < R.size(); i++){
+ for (unsigned int i = 0; i < R.size(); i++)
+ {
std::map<Dof, int>::iterator itR = unknown.find(R[i]);
if (itR != unknown.end()) NR[i] = itR->second;
else NR[i] = -1;
}
+ for (unsigned int i = 0; i < R.size(); i++)
+ {
+ if (NR[i] != -1)
+ {
+ _current->addToRightHandSide(NR[i], m(i));
+ }
+ }
+ }
- for (unsigned int i = 0; i < R.size(); i++){
- if (NR[i] != -1){
- for (unsigned int j = 0; j < R.size(); j++){
- if (NR[j] != -1){
- _current->addToMatrix(NR[i], NR[j], m(i, j));
- }
- else{
- typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(R[j]);
- if (itFixed != fixed.end()){
- _current->addToRightHandSide(NR[i], -m(i, j) * itFixed->second);
- }
- }
- }
+
+ inline void assemble(std::vector<Dof> &R, fullMatrix<dataMat> &m)
+ {
+ if (!_current->isAllocated()) _current->allocate(unknown.size());
+ std::vector<int> NR(R.size());
+ for (unsigned int i = 0; i < R.size(); i++)
+ {
+ std::map<Dof, int>::iterator itR = unknown.find(R[i]);
+ if (itR != unknown.end()) NR[i] = itR->second;
+ else NR[i] = -1;
+ }
+ for (unsigned int i = 0; i < R.size(); i++)
+ {
+ if (NR[i] != -1)
+ {
+ for (unsigned int j = 0; j < R.size(); j++)
+ {
+ if (NR[j] != -1)
+ {
+ _current->addToMatrix(NR[i], NR[j], m(i, j));
+ }
+ else
+ {
+ typename std::map<Dof, dataVec>::iterator itFixed = fixed.find(R[j]);
+ if (itFixed != fixed.end()){
+ _current->addToRightHandSide(NR[i], -m(i, j) * itFixed->second);
+ }
+ }
+ }
}
}
}
+
inline void assemble(int entR, int typeR, int entC, int typeC, const dataMat &value)
{
assemble(Dof(entR, typeR), Dof(entC, typeC), value);
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 35a5c62715..21b20c8590 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -313,13 +313,15 @@ void elasticitySolver::solve()
El.addToMatrix(*pAssembler, *elasticFields[i].g, *elasticFields[j].g);
}
}
-
+/*
printf("-- done assembling!\n");
- // for (int i=0;i<pAssembler->sizeOfR();i++){
- // printf("%g ",lsys->getFromRightHandSide(i));
- // }
- // printf("\n");
-
+ for (int i=0;i<pAssembler->sizeOfR();i++)
+ {
+ if (fabs(lsys->getFromRightHandSide(i))>0.000001)
+ printf("%g ",lsys->getFromRightHandSide(i));
+ }
+ printf("\n");
+*/
// solving
lsys->systemSolve();
printf("-- done solving!\n");
@@ -404,60 +406,109 @@ void MyelasticitySolver::solve()
}
}
- // line forces
- for (std::map<int, SVector3 >::iterator it = lineForces.begin();
- it != lineForces.end(); ++it){
- DummyfemTerm El(pModel);
+ VectorLagrangeFunctionSpace P123(_tag);//,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+
+ //line forces
+ for (std::map<int, SVector3 >::iterator it =lineForces.begin();
+ it != lineForces.end(); ++it)
+ {
int iEdge = it->first;
SVector3 f = it->second;
- El.neumannNodalBC(iEdge, 1, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
- El.neumannNodalBC(iEdge, 1, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
- El.neumannNodalBC(iEdge, 1, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
+ LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,f);
printf("-- Force on edge %3d : %8.5f %8.5f %8.5f\n", iEdge, f.x(), f.y(), f.z());
+ int dim=1;
+ std::map<int, std::vector<GEntity*> > groups[4];
+ pModel->getPhysicalGroups(groups);
+ fullVector<double> localVector;
+ std::vector<Dof> R;
+ std::map<int, std::vector<GEntity*> >::iterator itg = groups[dim].find(iEdge);
+ if (itg == groups[dim].end()) {printf(" Nonexistent edge\n");break;}
+ for (unsigned int i = 0; i < itg->second.size(); ++i)
+ {
+ GEntity *ge = itg->second[i];
+ for (unsigned int j = 0; j < ge->getNumMeshElements(); j++)
+ {
+ MElement *e = ge->getMeshElement(j);
+ int integrationOrder = 2 * e->getPolynomialOrder();
+ int npts;
+ IntPt *GP;
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ R.clear();
+ P123.getKeys(e,R);
+ Lterm.get(e,npts,GP,localVector);
+ pAssembler->assemble(R, localVector);
+ }
+ }
}
- // face forces
+
+//face forces
for (std::map<int, SVector3 >::iterator it = faceForces.begin();
- it != faceForces.end(); ++it){
- DummyfemTerm El(pModel);
+ it != faceForces.end(); ++it)
+ {
int iFace = it->first;
SVector3 f = it->second;
- El.neumannNodalBC(iFace, 2, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
- El.neumannNodalBC(iFace, 2, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
- El.neumannNodalBC(iFace, 2, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
+ LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,f);
printf("-- Force on face %3d : %8.5f %8.5f %8.5f\n", iFace, f.x(), f.y(), f.z());
+ int dim=2;
+ std::map<int, std::vector<GEntity*> > groups[4];
+ pModel->getPhysicalGroups(groups);
+ fullVector<double> localVector;
+ std::vector<Dof> R;
+ std::map<int, std::vector<GEntity*> >::iterator itg = groups[dim].find(iFace);
+ if (itg == groups[dim].end()) {printf(" Nonexistent face\n");break;}
+ for (unsigned int i = 0; i < itg->second.size(); ++i)
+ {
+ GEntity *ge = itg->second[i];
+ for (unsigned int j = 0; j < ge->getNumMeshElements(); j++)
+ {
+ MElement *e = ge->getMeshElement(j);
+ int integrationOrder = 2 * e->getPolynomialOrder();
+ int npts;
+ IntPt *GP;
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ R.clear();
+ P123.getKeys(e,R);
+ Lterm.get(e,npts,GP,localVector);
+ pAssembler->assemble(R, localVector);
+ }
+ }
}
- // volume forces
+
+//volume forces
for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
- it != volumeForces.end(); ++it){
- DummyfemTerm El(pModel);
+ it != volumeForces.end(); ++it)
+ {
int iVolume = it->first;
SVector3 f = it->second;
-// El.setVector(f);
+ LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,f);
printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n", iVolume, f.x(), f.y(), f.z());
- std::vector<GEntity*> ent = groups[_dim][iVolume];
- for (unsigned int i = 0; i < ent.size(); i++){
- // to do
- // El.addToRightHandSide(*pAssembler, ent[i]);
- }
- }
-
-/*
- // assembling the stifness matrix
- for (unsigned int i = 0; i < elasticFields.size(); i++){
- SElement::setShapeEnrichement(elasticFields[i]._enrichment);
- for (unsigned int j = 0; j < elasticFields.size(); j++){
- elasticityTerm El(0, elasticFields[i]._E, elasticFields[i]._nu,
- elasticFields[i]._tag, elasticFields[j]._tag);
- SElement::setTestEnrichement(elasticFields[j]._enrichment);
- // printf("%d %d\n",elasticFields[i]._tag,elasticFields[j]._tag);
- El.addToMatrix(*pAssembler, *elasticFields[i].g, *elasticFields[j].g);
+ int dim=3;
+ std::map<int, std::vector<GEntity*> > groups[4];
+ pModel->getPhysicalGroups(groups);
+ fullVector<double> localVector;
+ std::vector<Dof> R;
+ std::map<int, std::vector<GEntity*> >::iterator itg = groups[dim].find(iVolume);
+ if (itg == groups[dim].end()) {printf(" Nonexistent volume\n");break;}
+ for (unsigned int i = 0; i < itg->second.size(); ++i)
+ {
+ GEntity *ge = itg->second[i];
+ for (unsigned int j = 0; j < ge->getNumMeshElements(); j++)
+ {
+ MElement *e = ge->getMeshElement(j);
+ int integrationOrder = 2 * e->getPolynomialOrder();
+ int npts;
+ IntPt *GP;
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ R.clear();
+ P123.getKeys(e,R);
+ Lterm.get(e,npts,GP,localVector);
+ pAssembler->assemble(R, localVector);
+ }
}
}
-*/
- VectorLagrangeFunctionSpace P123(_tag);//,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
for (unsigned int i = 0; i < elasticFields.size(); i++)
{
@@ -477,6 +528,31 @@ void MyelasticitySolver::solve()
pAssembler->assemble(R, localMatrix);
}
}
+
+
+/*
+
+ for (std::map<int, SVector3 >::iterator it = volumeForces.begin();it != volumeForces.end(); ++it)
+ {
+ int iVolume = it->first;
+ SVector3 f = it->second;
+ printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n", iVolume, f.x(), f.y(), f.z());
+ std::vector<GEntity*> ent = groups[_dim][iVolume];
+ for (unsigned int i = 0; i < ent.size(); i++)
+ {
+ // to do
+ // El.addToRightHandSide(*pAssembler, ent[i]);
+ }
+ }
+*/
+/*
+ for (int i=0;i<pAssembler->sizeOfR();i++)
+ {
+ if (fabs(lsys->getFromRightHandSide(i))>0.000001)
+ printf("%g ",lsys->getFromRightHandSide(i));
+ }
+ printf("\n");
+*/
printf("-- done assembling!\n");
lsys->systemSolve();
printf("-- done solving!\n");
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 9007b27ae1..ca26866c02 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -106,7 +106,7 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
{
int ndofs= ele->getNumVertices();
int curpos=vals.size();
- vals.reserve(curpos+ndofs);
+ vals.resize(curpos+ndofs);
ele->getShapeFunctions(u, v, w, &(vals[curpos]));
};
virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
diff --git a/Solver/terms.h b/Solver/terms.h
index df76f2987d..b229a2f555 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -82,8 +82,11 @@ class LinearTermBase
template<class S1> class LinearTerm : public LinearTermBase
{
+ protected :
+ S1& space1;
public :
- virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m) =0;
+ LinearTerm(S1& space1_) : space1(space1_) {}
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m) =0;
};
class ScalarTermBase
@@ -168,6 +171,9 @@ template<class S1,class S2> class ElasticTerm : public BilinearTerm<S1,S2> // no
+
+
+
template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric
{
protected :
@@ -224,4 +230,34 @@ template<class S1> class ElasticTerm<S1,S1> : public BilinearTerm<S1,S1> // symm
}; // class
+inline double dot(const double &a, const double &b)
+{ return a*b; }
+
+template<class S1> class LoadTerm : public LinearTerm<S1>
+{
+ typename S1::ValType Load;
+ public :
+ LoadTerm(S1& space1_,typename S1::ValType Load_) :LinearTerm<S1>(space1_),Load(Load_) {};
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &m)
+ {
+ double nbFF=LinearTerm<S1>::space1.getNumKeys(ele);
+ double jac[3][3];
+ m.resize(nbFF);
+ m.scale(0.);
+ for (int i = 0; i < npts; i++)
+ {
+ const double u = GP[i].pt[0];const double v = GP[i].pt[1];const double w = GP[i].pt[2];
+ const double weight = GP[i].weight;const double detJ = ele->getJacobian(u, v, w, jac);
+ std::vector<typename S1::ValType> Vals;
+ LinearTerm<S1>::space1.f(ele,u, v, w, Vals);
+ for (int j = 0; j < nbFF ; ++j)
+ {
+ m(j)+=dot(Vals[j],Load)*weight*detJ;
+ }
+ }
+ }
+};
+
+
+
#endif// _TERMS_H_
--
GitLab