diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index 91bd3cf76bf3d97a4e9de4f0a41cf53f6f1f121c..c67eb0e072416cdab3e2347fc8bd153b10507355 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -74,10 +74,10 @@ class DofAffineConstraint{
// fullVecor<double>, ...
template <class T>
class dofManager{
- private:
+ public:
typedef typename dofTraits<T>::VecType dataVec;
typedef typename dofTraits<T>::MatType dataMat;
-
+ private:
// general affine constraint on sub-blocks, treated by adding
// equations:
// dataMat * DofVec = \sum_i dataMat_i * DofVec_i + dataVec
@@ -102,6 +102,10 @@ class dofManager{
public:
dofManager(linearSystem<dataMat> *l) : _current(l) { _linearSystems["A"] = l; }
+ inline void fixDof(Dof key, const dataVec &value)
+ {
+ fixed[key] = value;
+ }
inline void fixDof(long int ent, int type, const dataVec &value)
{
fixed[Dof(ent, type)] = value;
@@ -110,9 +114,8 @@ class dofManager{
{
fixDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField), value);
}
- inline void numberDof(long int ent, int type)
+ inline void numberDof(Dof key)
{
- Dof key(ent, type);
if(fixed.find(key) != fixed.end()) return;
// if (constraints.find(key) != constraints.end()) return;
std::map<Dof, int> :: iterator it = unknown.find(key);
@@ -121,6 +124,10 @@ class dofManager{
unknown[key] = size;
}
}
+ inline void numberDof(long int ent, int type)
+ {
+ numberDof(Dof(ent, type));
+ }
inline void numberVertex(MVertex*v, int iComp, int iField)
{
numberDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 1e25c25650253d07815f1de461393434ee8c5a64..38e33dacd1de6e2d48cf6c9bd0c35bc8ddb42133 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -16,6 +16,7 @@
#include "functionSpace.h"
#include "terms.h"
#include "solverAlgorithms.h"
+#include "quadratureRules.h"
#if defined(HAVE_POST)
#include "PView.h"
@@ -136,19 +137,31 @@ void elasticitySolver::readInputFile(const std::string &fn)
double val;
int node, comp;
if(fscanf(f, "%d %d %lf", &node, &comp, &val) != 3) return;
- nodalDisplacements[ std::make_pair(node, comp) ] = val;
+ nodalDisplacements[ std::make_pair(node, comp) ] = val;
}
else if (!strcmp(what, "EdgeDisplacement")){
double val;
int edge, comp;
if(fscanf(f, "%d %d %lf", &edge, &comp, &val) != 3) return;
- edgeDisplacements[ std::make_pair(edge, comp) ] = val;
+ edgeDisplacements[ std::make_pair(edge, comp) ] = val;
+ dirichletBC diri;
+ diri.g = new groupOfElements (1, edge);
+ diri._f= simpleFunction<double>(val);
+ diri._comp=comp;
+ diri._tag=edge;
+ edgeDiri.push_back(diri);
}
else if (!strcmp(what, "FaceDisplacement")){
double val;
int face, comp;
if(fscanf(f, "%d %d %lf", &face, &comp, &val) != 3) return;
- faceDisplacements[ std::make_pair(face, comp) ] = val;
+ faceDisplacements[ std::make_pair(face, comp) ] = val;
+ dirichletBC diri;
+ diri.g = new groupOfElements (2, face);
+ diri._f= simpleFunction<double>(val);
+ diri._comp=comp;
+ diri._tag=face;
+ faceDiri.push_back(diri);
}
else if (!strcmp(what, "NodalForce")){
double val1, val2, val3;
@@ -158,22 +171,37 @@ void elasticitySolver::readInputFile(const std::string &fn)
}
else if (!strcmp(what, "LineForce")){
double val1, val2, val3;
- int node;
- if(fscanf(f, "%d %lf %lf %lf", &node, &val1, &val2, &val3) != 4) return;
+ int edge;
+ if(fscanf(f, "%d %lf %lf %lf", &edge, &val1, &val2, &val3) != 4) return;
//printf("%d %lf %lf %lf\n", node, val1, val2, val3);
- lineForces[node] = SVector3(val1, val2, val3);
+ lineForces[edge] = SVector3(val1, val2, val3);
+ neumannBC neu;
+ neu.g = new groupOfElements (1, edge);
+ neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+ neu._tag=edge;
+ edgeNeu.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;
- faceForces[face] = SVector3(val1, val2, val3);
+ faceForces[face] = SVector3(val1, val2, val3);
+ neumannBC neu;
+ neu.g = new groupOfElements (2, face);
+ neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+ neu._tag=face;
+ edgeNeu.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;
- volumeForces[volume] = SVector3(val1, val2, val3);
+ volumeForces[volume] = SVector3(val1, val2, val3);
+ neumannBC neu;
+ neu.g = new groupOfElements (3, volume);
+ neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+ neu._tag=volume;
+ edgeNeu.push_back(neu);
}
else if (!strcmp(what, "MeshFile")){
char name[245];
@@ -349,6 +377,7 @@ void MyelasticitySolver::solve()
// 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
+ VectorLagrangeFunctionSpace P123(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
for (std::map<std::pair<int, int>, double>::iterator it = nodalDisplacements.begin();
it != nodalDisplacements.end(); ++it){
@@ -358,36 +387,24 @@ void MyelasticitySolver::solve()
it->first.first, it->first.second, it->second);
}
- for (std::map<std::pair<int, int>, double>::iterator it = edgeDisplacements.begin();
- it != edgeDisplacements.end(); ++it){
- DummyfemTerm El(pModel);
- El.dirichletNodalBC(it->first.first, 1, it->first.second, _tag,
- simpleFunction<double>(it->second), *pAssembler);
- printf("-- Fixing edge %3d comp %1d to %8.5f\n",
- it->first.first, it->first.second, it->second);
+
+ for (unsigned int i = 0; i < edgeDiri.size(); i++)
+ {
+ FixNodalDofs(P123,edgeDiri[i].g->begin(),edgeDiri[i].g->end(),*pAssembler,edgeDiri[i]._f,FilterDofComponent(edgeDiri[i]._comp));
+ printf("-- Fixing edge %3d comp %3d \n", edgeDiri[i]._tag, edgeDiri[i]._comp);
}
- for (std::map<std::pair<int, int>, double>::iterator it = faceDisplacements.begin();
- it != faceDisplacements.end(); ++it){
- DummyfemTerm El(pModel);
- El.dirichletNodalBC(it->first.first, 2, it->first.second, _tag,
- simpleFunction<double>(it->second), *pAssembler);
- printf("-- Fixing face %3d comp %1d to %8.5f\n",
- it->first.first, it->first.second, it->second);
+ for (unsigned int i = 0; i < faceDiri.size(); i++)
+ {
+ FixNodalDofs(P123,faceDiri[i].g->begin(),faceDiri[i].g->end(),*pAssembler,faceDiri[i]._f,FilterDofComponent(faceDiri[i]._comp));
+ printf("-- Fixing face %3d comp %3d \n", faceDiri[i]._tag, faceDiri[i]._comp);
}
// we number the nodes : when a node is numbered, it cannot be numbered
// again with another number. so we loop over elements
- for (unsigned int i = 0; i < elasticFields.size(); ++i){
- groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin();
- for ( ; it != elasticFields[i].g->end(); ++it){
- SElement se(*it);
- for (int j = 0; j < se.getNumNodalShapeFunctions(); ++j){
- pAssembler->numberVertex(se.getVertex(j), 0, elasticFields[i]._tag);
- pAssembler->numberVertex(se.getVertex(j), 1, elasticFields[i]._tag);
- pAssembler->numberVertex(se.getVertex(j), 2, elasticFields[i]._tag);
- }
- }
+ for (unsigned int i = 0; i < elasticFields.size(); ++i)
+ {
+ NumberDofs(P123, elasticFields[i].g->begin(), elasticFields[i].g->end(),*pAssembler);
}
// Now we start the assembly process
@@ -407,91 +424,37 @@ void MyelasticitySolver::solve()
}
}
- VectorLagrangeFunctionSpace P123(_tag);//,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
- //line forces
- for (std::map<int, SVector3 >::iterator it =lineForces.begin();
- it != lineForces.end(); ++it)
+ GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
+
+ for (unsigned int i = 0; i < edgeNeu.size(); i++)
{
- int iEdge = it->first;
- SVector3 f = it->second;
- 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*> >::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);
- Assemble(Lterm,P123,e,*pAssembler);
- }
- }
+ LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,edgeNeu[i]._f);
+ Assemble(Lterm,P123,edgeNeu[i].g->begin(),edgeNeu[i].g->end(),Integ_Boundary,*pAssembler);
+ printf("-- Force on edge %3d \n", edgeNeu[i]._tag);
}
-
-//face forces
- for (std::map<int, SVector3 >::iterator it = faceForces.begin();
- it != faceForces.end(); ++it)
+ for (unsigned int i = 0; i < faceNeu.size(); i++)
{
- int iFace = it->first;
- SVector3 f = it->second;
- 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*> >::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);
- Assemble(Lterm,P123,e,*pAssembler);
- }
- }
+ LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,faceNeu[i]._f);
+ Assemble(Lterm,P123,faceNeu[i].g->begin(),faceNeu[i].g->end(),Integ_Boundary,*pAssembler);
+ printf("-- Force on face %3d \n", faceNeu[i]._tag);
}
-
-//volume forces
- for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
- it != volumeForces.end(); ++it)
+ for (unsigned int i = 0; i < volumeNeu.size(); i++)
{
- int iVolume = it->first;
- SVector3 f = it->second;
- LoadTerm<VectorLagrangeFunctionSpace> Lterm(P123,f);
- printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n", iVolume, f.x(), f.y(), f.z());
- int dim=3;
- 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);
- Assemble(Lterm,P123,e,*pAssembler);
- }
- }
+ LoadTermSF<VectorLagrangeFunctionSpace> Lterm(P123,volumeNeu[i]._f);
+ Assemble(Lterm,P123,volumeNeu[i].g->begin(),volumeNeu[i].g->end(),Integ_Boundary,*pAssembler);
+ printf("-- Force on volume %3d \n", volumeNeu[i]._tag);
}
-
+ GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); i++)
{
ElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
- Assemble(Eterm,P123,elasticFields[i].g->begin(),elasticFields[i].g->end(),*pAssembler);
+ Assemble(Eterm,P123,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
}
-/*
- 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/elasticitySolver.h b/Solver/elasticitySolver.h
index c3bb8026d52f6f675a31becf2ae291e3018c9929..5a97d9cbb45adba80302732b4d193b399db72ab9 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -21,7 +21,22 @@ struct elasticField {
groupOfElements *g; // support for this field
simpleFunction<double> *_enrichment; // XFEM enrichment
double _E, _nu; // specific elastic datas (should be somewhere else)
- elasticField () : g(0), _enrichment(0){}
+ elasticField () : g(0), _enrichment(0),_tag(0){}
+};
+
+struct dirichletBC {
+ int _tag; // tag for the dofManager
+ int _comp; // component
+ groupOfElements *g; // support for this BC
+ simpleFunction<double> _f;
+ dirichletBC () : g(0),_comp(0),_tag(0){}
+};
+
+struct neumannBC {
+ int _tag; // tag for the dofManager
+ groupOfElements *g; // support for this BC
+ simpleFunction<SVector3> _f;
+ neumannBC () : g(0),_tag(0),_f(SVector3(0,0,0)){}
};
// an elastic solver ...
@@ -36,16 +51,22 @@ class elasticitySolver{
std::map<int, SVector3> nodalForces;
// imposed line forces
std::map<int, SVector3> lineForces;
+ std::vector<neumannBC> edgeNeu;
// imposed face forces
std::map<int, SVector3> faceForces;
+ std::vector<neumannBC> faceNeu;
// imposed volume forces
std::map<int, SVector3> volumeForces;
+ std::vector<neumannBC> volumeNeu;
// imposed nodal displacements
std::map<std::pair<int,int>, double> nodalDisplacements;
// imposed edge displacements
std::map<std::pair<int,int>, double> edgeDisplacements;
+ std::vector<dirichletBC> edgeDiri;
// imposed face displacements
std::map<std::pair<int,int>, double> faceDisplacements;
+ std::vector<dirichletBC> faceDiri;
+// std::vector<dirichletBC> faceDiri;
public:
elasticitySolver(int tag) : _tag(tag) {}
void addNodalForces (int iNode, const SVector3 &f)
diff --git a/Solver/quadratureRules.h b/Solver/quadratureRules.h
new file mode 100644
index 0000000000000000000000000000000000000000..4ae9438c0edeba1d25f41b4acd54a3d9252107da
--- /dev/null
+++ b/Solver/quadratureRules.h
@@ -0,0 +1,69 @@
+//
+// C++ Interface: quadratureRules
+//
+// Description:
+//
+//
+// Author: <Eric Bechet>, (C) 2009
+//
+// Copyright: See COPYING file that comes with this distribution
+//
+//
+
+
+#ifndef _QUADRATURERULES_H_
+#define _QUADRATURERULES_H_
+
+#include "Gauss.h"
+#include "MElement.h"
+
+class QuadratureBase
+{
+ public :
+ virtual ~QuadratureBase(){}
+ virtual int getIntPoints(MElement *e,IntPt **GP) =0;
+};
+
+
+class GaussQuadrature : public QuadratureBase
+{
+ public :
+ enum IntegCases {Other,Val,Grad,ValVal,GradGrad};
+ private :
+ int order;
+ IntegCases info;
+ public :
+ GaussQuadrature(int order_=0):order(order_),info(Other) {}
+ GaussQuadrature(IntegCases info_):order(0),info(info_) {}
+ virtual ~GaussQuadrature(){}
+ int getIntPoints(MElement *e,IntPt **GP)
+ {
+ int integrationOrder;
+ int npts;
+ int geoorder=e->getPolynomialOrder();
+ switch(info)
+ {
+ case Other :
+ integrationOrder = order;
+ break;
+ case Val :
+ integrationOrder=geoorder+1;
+ break;
+ case Grad :
+ integrationOrder=geoorder;
+ break;
+ case ValVal :
+ integrationOrder=2*geoorder;
+ break;
+ case GradGrad :
+ integrationOrder=3*(geoorder-1);
+ break;
+ default : integrationOrder=1;
+ }
+ e->getIntegrationPoints(integrationOrder, &npts, GP);
+ return npts;
+ }
+};
+
+
+#endif //_QUADRATURERULES_H_
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index 2337d0afcb62de4402e816591bb6fdd4de371b8b..a86e28117ae0bb184fb6348477f58d02188b2fe1 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -15,71 +15,140 @@
#define _SOLVERALGORITHMS_H_
-
+#include "dofManager.h"
#include "terms.h"
+#include "quadratureRules.h"
+#include "MVertex.h"
+#include <iostream>
+template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,QuadratureBase &integrator,Assembler &assembler) // symmetric
+{
+ fullMatrix<typename Assembler::dataMat> localMatrix;
+ std::vector<Dof> R;
+ for (Iterator it = itbegin;it!=itend; ++it)
+ {
+ MElement *e = *it;
+ R.clear();
+ IntPt *GP;
+ int npts=integrator.getIntPoints(e,&GP);
+ term.get(e,npts,GP,localMatrix);
+ space.getKeys(e,R);
+ assembler.assemble(R, localMatrix);
+ }
+}
-
-template<class Iterator,class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space1,Iterator itbegin,Iterator itend,Assembler &assembler) // symmetric
+template<class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler) // symmetric
{
- fullMatrix<double> localMatrix;
- std::vector<Dof> R;
- for (Iterator it = itbegin;it!=itend; ++it)
- {
- MElement *e = *it;
- R.clear();
- int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ; // todo: to be given from outside
- int npts=0;
- IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- term.get(e,npts,GP,localMatrix);
- space1.getKeys(e,R);
- assembler.assemble(R, localMatrix);
- }
+ fullMatrix<typename Assembler::dataMat> localMatrix;
+ std::vector<Dof> R;
+ IntPt *GP;
+ int npts=integrator.getIntPoints(e,&GP);
+ term.get(e,npts,GP,localMatrix);
+ space.getKeys(e,R);
+ assembler.assemble(R, localMatrix);
}
-template<class Assembler> void Assemble(BilinearTermBase &term,FunctionSpaceBase &space1,MElement *e,Assembler &assembler) // symmetric
+template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,QuadratureBase &integrator,Assembler &assembler)
{
- fullMatrix<double> localMatrix;
- std::vector<Dof> R;
- int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ; // todo: to be given from outside
- int npts=0;
+ fullVector<typename Assembler::dataMat> localVector;
+ std::vector<Dof> R;
+ for (Iterator it = itbegin;it!=itend; ++it)
+ {
+ MElement *e = *it;
+ R.clear();
IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- term.get(e,npts,GP,localMatrix);
- space1.getKeys(e,R);
- assembler.assemble(R, localMatrix);
+ int npts=integrator.getIntPoints(e,&GP);
+ term.get(e,npts,GP,localVector);
+ space.getKeys(e,R);
+ assembler.assemble(R, localVector);
+ }
+}
+
+template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,MElement *e,QuadratureBase &integrator,Assembler &assembler)
+{
+ fullVector<typename Assembler::dataMat> localVector;
+ std::vector<Dof> R;
+ IntPt *GP;
+ int npts=integrator.getIntPoints(e,&GP);
+ term.get(e,npts,GP,localVector);
+ space.getKeys(e,R);
+ assembler.assemble(R, localVector);
+}
+
+
+template<class Assembler> void FixDofs(Assembler &assembler,std::vector<Dof> &dofs,std::vector<typename Assembler::dataVec> &vals)
+{
+ int nbff=dofs.size();
+ for (int i=0;i<nbff;++i)
+ {
+ assembler.fixDof(dofs[i],vals[i]);
+ }
}
-template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space1,Iterator itbegin,Iterator itend,Assembler &assembler)
+class FilterDof
+{
+ public:
+ virtual bool operator()(Dof key) {return true;}
+};
+
+class FilterDofComponent :public FilterDof
{
- fullVector<double> localVector;
+ int comp;
+ public :
+ FilterDofComponent(int comp_):comp(comp_) {}
+ virtual bool operator()(Dof key)
+ {
+ int type=key.getType();
+ int icomp,iphys;
+ Dof::getTwoIntsFromType(type, iphys, icomp);
+ if (icomp==comp) return true;
+ return false;
+ }
+};
+
+
+template<class Iterator,class Assembler> void FixNodalDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof filter)
+{
+ for (Iterator it=itbegin;it!=itend;++it)
+ {
+ MElement *e=*it;
+ std::vector<MVertex*> tabV;
+ int nv=e->getNumVertices();
std::vector<Dof> R;
- for (Iterator it = itbegin;it!=itend; ++it)
+ space.getKeys(e,R);
+ tabV.reserve(nv);
+ for (int i=0;i<nv;++i) tabV.push_back(e->getVertex(i));
+
+ for (std::vector<Dof>::iterator itd=R.begin();itd!=R.end();++itd)
{
- MElement *e = *it;
- R.clear();
- int integrationOrder = 2* (e->getPolynomialOrder()) ; // todo: to be given from outside
- int npts=0;
- IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- term.get(e,npts,GP,localVector);
- space1.getKeys(e,R);
- assembler.assemble(R, localVector);
+ Dof key=*itd;
+ if (filter(key))
+ {
+ for (int i=0;i<nv;++i)
+ {
+ if (tabV[i]->getNum()==key.getEntity())
+ {
+ assembler.fixDof(key, fct(tabV[i]->x(),tabV[i]->y(),tabV[i]->z()));
+ break;
+ }
+ }
+ }
}
+ }
}
-template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space1,MElement *e,Assembler &assembler)
+
+template<class Iterator,class Assembler> void NumberDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler)
{
- fullVector<double> localVector;
+ for (Iterator it=itbegin;it!=itend;++it)
+ {
+ MElement *e=*it;
std::vector<Dof> R;
- int integrationOrder = 2* (e->getPolynomialOrder()) ; // todo: to be given from outside
- int npts=0;
- IntPt *GP;
- e->getIntegrationPoints(integrationOrder, &npts, &GP);
- term.get(e,npts,GP,localVector);
- space1.getKeys(e,R);
- assembler.assemble(R, localVector);
+ space.getKeys(e,R);
+ int nbdofs=R.size();
+ for (int i=0;i<nbdofs;++i) assembler.numberDof(R[i]);
+ }
}
-#endif// _SOLVERALGORITHMS_H_
\ No newline at end of file
+
+#endif// _SOLVERALGORITHMS_H_
diff --git a/Solver/terms.h b/Solver/terms.h
index ab5ab4652fc9deb571a96948263fdebbee1cdaa0..6ac36ff7d6de9807e50d1af3e09177580162c411 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -19,6 +19,7 @@
#include <iterator>
#include "Numeric.h"
#include "functionSpace.h"
+#include "groupOfElements.h"
// evaluation of a field ???
@@ -259,5 +260,32 @@ template<class S1> class LoadTerm : public LinearTerm<S1>
}
};
+template<class S1> class LoadTermSF : public LinearTerm<S1>
+{
+ simpleFunction<typename S1::ValType> Load;
+ public :
+ LoadTermSF(S1& space1_,simpleFunction<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(u,v,w))*weight*detJ;
+ }
+ }
+ }
+};
+
+
+
#endif// _TERMS_H_