From 6146be39b0f51e00777712cddfc86512820618da Mon Sep 17 00:00:00 2001
From: Eric Bechet <eric.bechet@ulg.ac.be>
Date: Wed, 16 Dec 2009 17:38:09 +0000
Subject: [PATCH] weight watchers
---
Solver/dofManager.h | 22 ++
Solver/elasticitySolver.cpp | 431 ++++++++++++------------------------
Solver/elasticitySolver.h | 82 +++----
Solver/functionSpace.h | 105 +++++++--
Solver/solverAlgorithms.h | 79 +++++--
Solver/terms.h | 81 ++++++-
6 files changed, 415 insertions(+), 385 deletions(-)
diff --git a/Solver/dofManager.h b/Solver/dofManager.h
index c67eb0e072..40bc19e3ff 100644
--- a/Solver/dofManager.h
+++ b/Solver/dofManager.h
@@ -132,6 +132,28 @@ class dofManager{
{
numberDof(v->getNum(), Dof::createTypeWithTwoInts(iComp, iField));
}
+
+ inline void getDofValue(std::vector<Dof> &keys,std::vector<dataVec> &Vals)
+ {
+ int ndofs=keys.size();
+ Vals.reserve(Vals.size()+ndofs);
+ for (int i=0;i<ndofs;++i) Vals.push_back(getDofValue(keys[i]));
+ }
+
+ inline dataVec getDofValue(Dof key) const
+ {
+ {
+ typename std::map<Dof, dataVec>::const_iterator it = fixed.find(key);
+ if (it != fixed.end()) return it->second;
+ }
+ {
+ std::map<Dof, int>::const_iterator it = unknown.find(key);
+ if (it != unknown.end())
+ return _current->getFromSolution(it->second);
+ }
+ return dataVec(0.0);
+ }
+
inline dataVec getDofValue(int ent, int type) const
{
Dof key(ent, type);
diff --git a/Solver/elasticitySolver.cpp b/Solver/elasticitySolver.cpp
index 2d2599c0b5..e7e5575cae 100644
--- a/Solver/elasticitySolver.cpp
+++ b/Solver/elasticitySolver.cpp
@@ -5,14 +5,12 @@
#include <string.h>
#include "GmshConfig.h"
-#include "elasticityTerm.h"
-#include "MTriangle.h"
-#include "MTetrahedron.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"
@@ -21,93 +19,17 @@
#if defined(HAVE_POST)
#include "PView.h"
#include "PViewData.h"
-
-PView* elasticitySolver::buildDisplacementView (const std::string &postFileName){
- std::map<int, std::vector<double> > data;
- std::set<MVertex*> v;
- 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){
- v.insert(se.getVertex(j));
- }
- }
- }
-
- std::set<MVertex*>::iterator it = v.begin();
- for ( ; it != v.end(); ++it){
- std::vector<double> d;
- d.push_back(0); d.push_back(0); d.push_back(0);
- for (unsigned int i = 0; i < elasticFields.size(); ++i){
- const double E = (elasticFields[i]._enrichment) ?
- (*elasticFields[i]._enrichment)((*it)->x(), (*it)->y(), (*it)->z()) : 1.0;
- // printf("%g %d \n",pAssembler->getDofValue(*it,0,elasticFields[i]._tag),elasticFields[i]._tag);
- d[0] += pAssembler->getDofValue(*it, 0, elasticFields[i]._tag) * E;
- d[1] += pAssembler->getDofValue(*it, 1, elasticFields[i]._tag) * E;
- d[2] += pAssembler->getDofValue(*it, 2, elasticFields[i]._tag) * E;
- }
- data[(*it)->getNum()] = d;
- }
-
- PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
- return pv;
-}
-
-#else
-PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
-{
- Msg::Error("Post-pro module not available");
- return 0;
-}
#endif
-
-static double vonMises(dofManager<double> *a, MElement *e,
- double u, double v, double w,
- double E, double nu, int _tag)
-{
- double valx[256];
- double valy[256];
- double valz[256];
- for (int k = 0; k < e->getNumVertices(); k++){
- valx[k] = a->getDofValue(e->getVertex(k), 0, _tag);
- valy[k] = a->getDofValue(e->getVertex(k), 1, _tag);
- valz[k] = a->getDofValue(e->getVertex(k), 2, _tag);
- }
- double gradux[3];
- double graduy[3];
- double graduz[3];
- e->interpolateGrad(valx, u, v, w, gradux);
- e->interpolateGrad(valy, u, v, w, graduy);
- e->interpolateGrad(valz, u, v, w, graduz);
-
- double eps[6] = {gradux[0], graduy[1], graduz[2],
- 0.5 * (gradux[1] + graduy[0]),
- 0.5 * (gradux[2] + graduz[0]),
- 0.5 * (graduy[2] + graduz[1])};
- double A = E / (1. + nu);
- double B = A * (nu / (1. - 2 * nu));
- double trace = eps[0] + eps[1] + eps[2] ;
- double sxx = A * eps[0] + B * trace;
- double syy = A * eps[1] + B * trace;
- double szz = A * eps[2] + B * trace;
- double sxy = A * eps[3];
- double sxz = A * eps[4];
- double syz = A * eps[5];
-
- double s[9] = {sxx, sxy, sxz,
- sxy, syy, syz,
- sxz, syz, szz};
-
- return ComputeVonMises(s);
-}
-
void elasticitySolver::setMesh(const std::string &meshFileName)
{
pModel = new GModel();
pModel->readMSH(meshFileName.c_str());
_dim = pModel->getNumRegions() ? 3 : 2;
+ if (LagSpace) delete LagSpace;
+ if (_dim==3) LagSpace=new VectorLagrangeFunctionSpace(_tag);
+ if (_dim==2) LagSpace=new VectorLagrangeFunctionSpace(_tag,VectorLagrangeFunctionSpace::VECTOR_X,VectorLagrangeFunctionSpace::VECTOR_Y);
+
}
void elasticitySolver::readInputFile(const std::string &fn)
@@ -116,13 +38,11 @@ void elasticitySolver::readInputFile(const std::string &fn)
char what[256];
while(!feof(f)){
if(fscanf(f, "%s", what) != 1) return;
- // printf("%s\n", what);
if (!strcmp(what, "ElasticDomain")){
elasticField field;
int physical;
if(fscanf(f, "%d %lf %lf", &physical, &field._E, &field._nu) != 3) return;
field._tag = _tag;
- // printf("E %g nu %g\n",field._E,field._nu);
field.g = new groupOfElements (_dim, physical);
elasticFields.push_back(field);
}
@@ -133,75 +53,85 @@ void elasticitySolver::readInputFile(const std::string &fn)
// assign a tag
// elasticFields.push_back(field);
}
- else if (!strcmp(what, "NodalDisplacement")){
+ else if (!strcmp(what, "NodeDisplacement")){
double val;
int node, comp;
if(fscanf(f, "%d %d %lf", &node, &comp, &val) != 3) return;
- nodalDisplacements[ std::make_pair(node, comp) ] = val;
+ dirichletBC diri;
+ diri.g = new groupOfElements (0, node);
+ diri._f= simpleFunction<double>(val);
+ diri._comp=comp;
+ diri._tag=node;
+ diri.onWhat=BoundaryCondition::ON_VERTEX;
+ allDirichlet.push_back(diri);
}
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;
dirichletBC diri;
diri.g = new groupOfElements (1, edge);
diri._f= simpleFunction<double>(val);
diri._comp=comp;
diri._tag=edge;
- edgeDiri.push_back(diri);
+ diri.onWhat=BoundaryCondition::ON_EDGE;
+ allDirichlet.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;
dirichletBC diri;
diri.g = new groupOfElements (2, face);
diri._f= simpleFunction<double>(val);
diri._comp=comp;
diri._tag=face;
- faceDiri.push_back(diri);
+ diri.onWhat=BoundaryCondition::ON_FACE;
+ allDirichlet.push_back(diri);
}
- else if (!strcmp(what, "NodalForce")){
+ else if (!strcmp(what, "NodeForce")){
double val1, val2, val3;
int node;
if(fscanf(f, "%d %lf %lf %lf", &node, &val1, &val2, &val3) != 4) return;
- nodalForces[node] = SVector3(val1, val2, val3);
+ neumannBC neu;
+ neu.g = new groupOfElements (0, node);
+ neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
+ neu._tag=node;
+ neu.onWhat=BoundaryCondition::ON_VERTEX;
+ allNeumann.push_back(neu);
}
- else if (!strcmp(what, "LineForce")){
+ else if (!strcmp(what, "EdgeForce")){
double val1, val2, val3;
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[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);
+ neu.onWhat=BoundaryCondition::ON_EDGE;
+ allNeumann.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);
neumannBC neu;
neu.g = new groupOfElements (2, face);
neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
neu._tag=face;
- edgeNeu.push_back(neu);
+ neu.onWhat=BoundaryCondition::ON_FACE;
+ allNeumann.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);
neumannBC neu;
neu.g = new groupOfElements (3, volume);
neu._f= simpleFunction<SVector3>(SVector3(val1, val2, val3));
neu._tag=volume;
- edgeNeu.push_back(neu);
+ neu.onWhat=BoundaryCondition::ON_VOLUME;
+ allNeumann.push_back(neu);
}
else if (!strcmp(what, "MeshFile")){
char name[245];
@@ -210,7 +140,7 @@ void elasticitySolver::readInputFile(const std::string &fn)
}
else {
Msg::Error("Invalid input : %s", what);
- return;
+// return;
}
}
fclose(f);
@@ -226,236 +156,155 @@ void elasticitySolver::solve()
linearSystemGmm<double> *lsys = new linearSystemGmm<double>;
lsys->setNoisy(2);
#endif
+ if (pAssembler) delete pAssembler;
pAssembler = new dofManager<double>(lsys);
- std::map<int, std::vector<GEntity*> > groups[4];
- pModel->getPhysicalGroups(groups);
-
// 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 (std::map<std::pair<int, int>, double>::iterator it = nodalDisplacements.begin();
- it != nodalDisplacements.end(); ++it){
- MVertex *v = pModel->getMeshVertexByTag(it->first.first);
- pAssembler->fixVertex(v, it->first.second, _tag, it->second);
- printf("-- Fixing node %3d comp %1d to %8.5f\n",
- it->first.first, it->first.second, it->second);
- }
-
- for (std::map<std::pair<int, int>, double>::iterator it = edgeDisplacements.begin();
- it != edgeDisplacements.end(); ++it){
- elasticityTerm El(pModel, 1, 1, _tag);
- 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 (std::map<std::pair<int, int>, double>::iterator it = faceDisplacements.begin();
- it != faceDisplacements.end(); ++it){
- elasticityTerm El(pModel, 1, 1, _tag);
- 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 < 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 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);
- }
- }
+ // 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
- // Nodes at geometrical vertices
- for (std::map<int, SVector3 >::iterator it = nodalForces.begin();
- it != nodalForces.end(); ++it){
- int iVertex = it->first;
- SVector3 f = it->second;
- std::vector<GEntity*> ent = groups[0][iVertex];
- for (unsigned int i = 0; i < ent.size(); i++){
- pAssembler->assemble(ent[i]->mesh_vertices[0], 0, _tag, f.x());
- pAssembler->assemble(ent[i]->mesh_vertices[0], 1, _tag, f.y());
- pAssembler->assemble(ent[i]->mesh_vertices[0], 2, _tag, f.z());
- printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n",
- ent[i]->mesh_vertices[0]->getNum(), iVertex, f.x(), f.y(), f.z());
- }
- }
- // line forces
- for (std::map<int, SVector3 >::iterator it = lineForces.begin();
- it != lineForces.end(); ++it){
- elasticityTerm El(pModel, 1, 1, _tag);
- 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);
- printf("-- Force on edge %3d : %8.5f %8.5f %8.5f\n", iEdge, f.x(), f.y(), f.z());
- }
+ GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
- // face forces
- for (std::map<int, SVector3 >::iterator it = faceForces.begin();
- it != faceForces.end(); ++it){
- elasticityTerm El(pModel, 1, 1, _tag);
- 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);
- printf("-- Force on face %3d : %8.5f %8.5f %8.5f\n", iFace, f.x(), f.y(), f.z());
+ for (unsigned int i = 0; i < allNeumann.size(); i++)
+ {
+ LoadTerm<FunctionSpace<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);
}
- // volume forces
- for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
- it != volumeForces.end(); ++it){
- elasticityTerm El(pModel, 1, 1, _tag);
- int iVolume = it->first;
- SVector3 f = it->second;
- El.setVector(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]);
- }
- }
+// bulk material law
- // 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);
- }
+ GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
+ for (unsigned int i = 0; i < elasticFields.size(); i++)
+ {
+ IsotropicElasticTerm<FunctionSpace<SVector3> ,FunctionSpace<SVector3> > 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");
- 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");
}
+#if defined(HAVE_POST)
-void MyelasticitySolver::solve()
+static double vonMises(dofManager<double> *a, MElement *e,
+ double u, double v, double w,
+ double E, double nu, int _tag)
{
-#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
- pAssembler = new dofManager<double>(lsys);
-
- std::map<int, std::vector<GEntity*> > groups[4];
- pModel->getPhysicalGroups(groups);
-
- // 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){
- MVertex *v = pModel->getMeshVertexByTag(it->first.first);
- pAssembler->fixVertex(v, it->first.second, _tag, it->second);
- printf("-- Fixing node %3d comp %1d to %8.5f\n",
- it->first.first, it->first.second, it->second);
+ double valx[256];
+ double valy[256];
+ double valz[256];
+ for (int k = 0; k < e->getNumVertices(); k++){
+ valx[k] = a->getDofValue(e->getVertex(k), 0, _tag);
+ valy[k] = a->getDofValue(e->getVertex(k), 1, _tag);
+ valz[k] = a->getDofValue(e->getVertex(k), 2, _tag);
}
+ double gradux[3];
+ double graduy[3];
+ double graduz[3];
+ e->interpolateGrad(valx, u, v, w, gradux);
+ e->interpolateGrad(valy, u, v, w, graduy);
+ e->interpolateGrad(valz, u, v, w, graduz);
+ double eps[6] = {gradux[0], graduy[1], graduz[2],
+ 0.5 * (gradux[1] + graduy[0]),
+ 0.5 * (gradux[2] + graduz[0]),
+ 0.5 * (graduy[2] + graduz[1])};
+ double A = E / (1. + nu);
+ double B = A * (nu / (1. - 2 * nu));
+ double trace = eps[0] + eps[1] + eps[2] ;
+ double sxx = A * eps[0] + B * trace;
+ double syy = A * eps[1] + B * trace;
+ double szz = A * eps[2] + B * trace;
+ double sxy = A * eps[3];
+ double sxz = A * eps[4];
+ double syz = A * eps[5];
- 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);
- }
+ double s[9] = {sxx, sxy, sxz,
+ sxy, syy, syz,
+ sxz, syz, szz};
- 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);
- }
+ return ComputeVonMises(s);
+}
- // we number the nodes : when a node is numbered, it cannot be numbered
- // again with another number. so we loop over elements
+PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
+{
+ std::set<MVertex*> v;
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
- // First build the force vector
- // Nodes at geometrical vertices
- for (std::map<int, SVector3 >::iterator it = nodalForces.begin();
- it != nodalForces.end(); ++it){
- int iVertex = it->first;
- SVector3 f = it->second;
- std::vector<GEntity*> ent = groups[0][iVertex];
- for (unsigned int i = 0; i < ent.size(); i++){
- pAssembler->assemble(ent[i]->mesh_vertices[0], 0, _tag, f.x());
- pAssembler->assemble(ent[i]->mesh_vertices[0], 1, _tag, f.y());
- pAssembler->assemble(ent[i]->mesh_vertices[0], 2, _tag, f.z());
- printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n",
- ent[i]->mesh_vertices[0]->getNum(), iVertex, f.x(), f.y(), f.z());
+ for (groupOfElements::elementContainer::const_iterator it = elasticFields[i].g->begin(); it != elasticFields[i].g->end(); ++it)
+ {
+ MElement *e=*it;
+ for (int j = 0; j < e->getNumVertices(); ++j) v.insert(e->getVertex(j));
}
}
-
-
- GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
-
- for (unsigned int i = 0; i < edgeNeu.size(); i++)
+ std::map<int, std::vector<double> > data;
+ for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
{
- LoadTerm<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);
+ SVector3 val(0,0,0);
+ for (unsigned int i = 0; i < elasticFields.size(); ++i)
+ {
+ std::vector<Dof> D;
+ std::vector<SVector3> SFVals;
+ std::vector<double> Vals;
+ LagSpace->getKeys(*it,D);
+ pAssembler->getDofValue(D,Vals);
+ LagSpace->f(*it,SFVals);
+ for (int i=0;i<D.size();++i)
+ val+=SFVals[i]*Vals[i];
+ }
+ std::vector<double> vec(3);vec[0]=val(0);vec[1]=val(1);vec[2]=val(2);
+ data[(*it)->getNum()]=vec;
}
+ PView *pv = new PView (postFileName, "NodeData", pModel, data, 0.0);
+ return pv;
+}
- for (unsigned int i = 0; i < faceNeu.size(); i++)
- {
- LoadTerm<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);
- }
+PView *elasticitySolver::buildVonMisesView(const std::string &postFileName)
+{
+ std::map<int, std::vector<double> > data;
+ PView *pv = new PView (postFileName, "ElementData", pModel, data, 0.0);
+ return pv;
+}
- for (unsigned int i = 0; i < volumeNeu.size(); i++)
- {
- LoadTerm<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++)
- {
- IsotropicElasticTerm<VectorLagrangeFunctionSpace,VectorLagrangeFunctionSpace> Eterm(P123,elasticFields[i]._E,elasticFields[i]._nu);
- Assemble(Eterm,P123,elasticFields[i].g->begin(),elasticFields[i].g->end(),Integ_Bulk,*pAssembler);
- }
+#else
+PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
+{
+ Msg::Error("Post-pro module not available");
+ return 0;
+}
- printf("-- done assembling!\n");
- lsys->systemSolve();
- printf("-- done solving!\n");
+PView* elasticitySolver::buildVonMisesView(const std::string &postFileName)
+{
+ Msg::Error("Post-pro module not available");
+ return 0;
}
+
+#endif
+
diff --git a/Solver/elasticitySolver.h b/Solver/elasticitySolver.h
index 5a97d9cbb4..d9b82a8fe7 100644
--- a/Solver/elasticitySolver.h
+++ b/Solver/elasticitySolver.h
@@ -11,6 +11,7 @@
#include "SVector3.h"
#include "dofManager.h"
#include "simpleFunction.h"
+#include "functionSpace.h"
class GModel;
class PView;
@@ -24,68 +25,56 @@ struct elasticField {
elasticField () : g(0), _enrichment(0),_tag(0){}
};
-struct dirichletBC {
+struct BoundaryCondition
+{
+ enum location{UNDEF,ON_VERTEX,ON_EDGE,ON_FACE,ON_VOLUME};
+ location onWhat; // on vertices or elements
int _tag; // tag for the dofManager
- int _comp; // component
groupOfElements *g; // support for this BC
+ BoundaryCondition() : g(0),_tag(0),onWhat(UNDEF) {};
+};
+
+struct dirichletBC : public BoundaryCondition
+{
+ int _comp; // component
simpleFunction<double> _f;
- dirichletBC () : g(0),_comp(0),_tag(0){}
+ dirichletBC () :BoundaryCondition(),_comp(0),_f(0){}
};
-struct neumannBC {
- int _tag; // tag for the dofManager
- groupOfElements *g; // support for this BC
+struct neumannBC : public BoundaryCondition
+{
simpleFunction<SVector3> _f;
- neumannBC () : g(0),_tag(0),_f(SVector3(0,0,0)){}
+ neumannBC () : BoundaryCondition(),_f(SVector3(0,0,0)){}
};
// an elastic solver ...
-class elasticitySolver{
+class elasticitySolver
+{
protected:
GModel *pModel;
int _dim, _tag;
dofManager<double> *pAssembler;
+ FunctionSpace<SVector3> *LagSpace;
+
// young modulus and poisson coefficient per physical
std::vector<elasticField> elasticFields;
- // imposed nodal forces
- 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;
+ // neumann BC
+ std::vector<neumannBC> allNeumann;
+ // dirichlet BC
+ std::vector<dirichletBC> allDirichlet;
+
public:
- elasticitySolver(int tag) : _tag(tag) {}
- void addNodalForces (int iNode, const SVector3 &f)
+ elasticitySolver(int tag) : _tag(tag),LagSpace(0),pAssembler(0) {}
+ virtual ~elasticitySolver()
{
- nodalForces[iNode] = f;
+ if (LagSpace) delete LagSpace;
+ if (pAssembler) delete pAssembler;
}
- void addNodalDisplacement(int iNode, int dir, double val)
- {
- nodalDisplacements[std::make_pair(iNode, dir)] = val;
- }
- // void addElasticConstants(double e, double nu, int physical)
- // {
- // elasticConstants[physical] = std::make_pair(e, nu);
- // }
+ void readInputFile(const std::string &meshFileName);
void setMesh(const std::string &meshFileName);
virtual void solve();
- void readInputFile(const std::string &meshFileName);
virtual PView *buildDisplacementView(const std::string &postFileName);
- // PView *buildVonMisesView(const std::string &postFileName);
+ virtual PView *buildVonMisesView(const std::string &postFileName);
// std::pair<PView *, PView*> buildErrorEstimateView
// (const std::string &errorFileName, double, int);
// std::pair<PView *, PView*> buildErrorEstimateView
@@ -93,15 +82,4 @@ class elasticitySolver{
};
-
-// another elastic solver ...
-class MyelasticitySolver : public elasticitySolver
-{
- public:
- MyelasticitySolver(int tag) : elasticitySolver(tag) {}
- virtual void solve();
-};
-
-
-
#endif
diff --git a/Solver/functionSpace.h b/Solver/functionSpace.h
index 105fea2127..659037c2f4 100644
--- a/Solver/functionSpace.h
+++ b/Solver/functionSpace.h
@@ -67,13 +67,15 @@ class FunctionSpaceBase
{
public:
virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
- virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
+ virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
+ virtual int getNumKeys(MVertex *ver)=0; // if one needs the number of dofs
+ virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)=0;
};
template<class T>
class FunctionSpace : public FunctionSpaceBase
{
- protected:
+ public:
typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType;
/* typedef typename TensorialTraits<T>::HessType HessType;
@@ -81,6 +83,7 @@ class FunctionSpace : public FunctionSpaceBase
typedef typename TensorialTraits<T>::CurlType CurlType;*/
public:
+ virtual int f(MVertex *ver, std::vector<ValType> &vals) {} // used for neumann BC
virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)=0;
virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)=0;
// virtual groupOfElements* getSupport()=0;// probablement inutile
@@ -89,26 +92,24 @@ class FunctionSpace : public FunctionSpaceBase
virtual int divf(MElement *ele, double u, double v, double w,std::vector<DivType> &divs);
virtual int curlf(MElement *ele, double u, double v, double w,std::vector<CurlType> &curls);*/
virtual int getNumKeys(MElement *ele)=0; // if one needs the number of dofs
- virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
+ virtual int getKeys(MElement *ele, std::vector<Dof> &keys)=0;
+ virtual int getNumKeys(MVertex *ver)=0; // if one needs the number of dofs
+ virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)=0;
};
class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
{
- protected:
+ 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::vector<basisFunction*> basisFunctions;
int _iField; // field number (used to build dof keys)
- Dof getLocalDof(MElement *ele, int i) const
- {
- return Dof(ele->getVertex(i)->getNum(), _iField);
- }
-
public:
ScalarLagrangeFunctionSpace(int i=0):_iField(i) {}
virtual int getId(void) const {return _iField;};
@@ -118,7 +119,8 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
int curpos=vals.size();
vals.resize(curpos+ndofs);
ele->getShapeFunctions(u, v, w, &(vals[curpos]));
- };
+ };
+ virtual int f(MVertex *ver, std::vector<ValType> &vals) {vals.push_back(1.0);} // used for neumann BC
virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads)
{
int ndofs= ele->getNumVertices();
@@ -143,8 +145,13 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
int ndofs= ele->getNumVertices();
keys.reserve(keys.size()+ndofs);
for (int i=0;i<ndofs;++i)
- keys.push_back(getLocalDof(ele,i));
+ getKeys(ele->getVertex(i),keys);
}
+ virtual int getNumKeys(MVertex *ver) { return 1;}
+ virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)
+ {
+ keys.push_back(Dof(ver->getNum(), _iField));
+ };
};
@@ -180,6 +187,20 @@ public :
}
virtual ~ScalarToAnyFunctionSpace() {delete ScalarFS;}
+ virtual int f(MVertex *ver, std::vector<ValType> &vals)
+ {
+ std::vector<double> valsd;
+ ScalarFS->f(ver,valsd);
+ int nbdofs=valsd.size();
+ int nbcomp=comp.size();
+ int curpos=vals.size();
+ vals.reserve(curpos+nbcomp*nbdofs);
+ for (int j=0;j<nbcomp;++j)
+ {
+ for (int i=0;i<nbdofs;++i) vals.push_back(multipliers[j]*valsd[i]);
+ }
+ } // used for neumann BC
+
virtual int f(MElement *ele, double u, double v, double w, std::vector<ValType> &vals)
{
std::vector<double> valsd;
@@ -236,14 +257,36 @@ public :
}
}
}
+
+ virtual int getNumKeys(MVertex *ver) { return ScalarFS->getNumKeys(ver)*comp.size();}
+ virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)
+ {
+ int nk=ScalarFS->getNumKeys(ver);
+ std::vector<Dof> bufk;
+ bufk.reserve(nk);
+ ScalarFS->getKeys(ver,bufk);
+ int nbdofs=bufk.size();
+ int nbcomp=comp.size();
+ int curpos=keys.size();
+ keys.reserve(curpos+nbcomp*nbdofs);
+ for (int j=0;j<nbcomp;++j)
+ {
+ for (int i=0;i<nk;++i)
+ {
+ int i1,i2;
+ Dof::getTwoIntsFromType(bufk[i].getType(), i1,i2);
+ keys.push_back(Dof(bufk[i].getEntity(),Dof::createTypeWithTwoInts(comp[j],i1)));
+ }
+ }
+ };
};
class VectorLagrangeFunctionSpace : public ScalarToAnyFunctionSpace<SVector3>
{
+ protected:
+ static const SVector3 BasisVectors[3];
public:
enum Along { VECTOR_X=0, VECTOR_Y=1, VECTOR_Z=2 };
- static const SVector3 BasisVectors[3];
-
typedef TensorialTraits<SVector3>::ValType ValType;
typedef TensorialTraits<SVector3>::GradType GradType;
typedef TensorialTraits<SVector3>::HessType HessType;
@@ -305,6 +348,12 @@ class CompositeFunctionSpace : public FunctionSpace<T>
delete (*it);
}
+ virtual int f(MVertex *ver, std::vector<ValType> &vals)
+ {
+ for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
+ (*it)->f(ver,vals);
+ }
+
virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals)
{
for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
@@ -330,11 +379,26 @@ class CompositeFunctionSpace : public FunctionSpace<T>
for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
(*it)->getKeys(ele,keys);
}
+
+ virtual int getNumKeys(MVertex *ver)
+ {
+ int ndofs=0;
+ for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
+ ndofs+=(*it)->getNumKeys(ver);
+ return ndofs;
+ }
+
+ virtual int getKeys(MVertex *ver, std::vector<Dof> &keys)
+ {
+ for (iterFS it=_spaces.begin(); it!=_spaces.end();++it)
+ (*it)->getKeys(ver,keys);
+ };
};
template<class T>
class xFemFunctionSpace : public FunctionSpace<T>
{
+ public:
typedef typename TensorialTraits<T>::ValType ValType;
typedef typename TensorialTraits<T>::GradType GradType;
typedef typename TensorialTraits<T>::HessType HessType;
@@ -344,10 +408,14 @@ class xFemFunctionSpace : public FunctionSpace<T>
FunctionSpace<T>* _spacebase;
// Function<double>* enrichment;
public:
+ virtual int f(MVertex *ver, std::vector<ValType> &vals);
virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
- int getNumKeys(MElement *ele);
- void getKeys(MElement *ele, std::vector<Dof> &keys);
+ virtual int getNumKeys(MElement *ele);
+ virtual void getKeys(MElement *ele, std::vector<Dof> &keys);
+ virtual int getNumKeys(MVertex *ver);
+ virtual int getKeys(MVertex *ver, std::vector<Dof> &keys);
+
};
@@ -364,10 +432,13 @@ class FilteredFunctionSpace : public FunctionSpace<T>
F &_filter;
public:
+ virtual int f(MVertex *ver, std::vector<ValType> &vals);
virtual int f(MElement *ele, double u, double v, double w,std::vector<ValType> &vals);
virtual int gradf(MElement *ele, double u, double v, double w,std::vector<GradType> &grads);
- int getNumKeys(MElement *ele);
- void getKeys(MElement *ele, std::vector<Dof> &keys);
+ virtual int getNumKeys(MElement *ele);
+ virtual void getKeys(MElement *ele, std::vector<Dof> &keys);
+ virtual int getNumKeys(MVertex *ver);
+ virtual int getKeys(MVertex *ver, std::vector<Dof> &keys);
};
diff --git a/Solver/solverAlgorithms.h b/Solver/solverAlgorithms.h
index 70c1abb848..0afc37ef51 100644
--- a/Solver/solverAlgorithms.h
+++ b/Solver/solverAlgorithms.h
@@ -21,6 +21,7 @@
#include "MVertex.h"
+
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;
@@ -64,6 +65,20 @@ template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,Func
}
}
+template<class Iterator,class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler)
+{
+ fullVector<typename Assembler::dataMat> localVector;
+ std::vector<Dof> R;
+ for (Iterator it = itbegin;it!=itend; ++it)
+ {
+ MVertex *v = *it;
+ R.clear();
+ term.get(v,localVector);
+ space.getKeys(v,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;
@@ -75,7 +90,6 @@ template<class Assembler> void Assemble(LinearTermBase &term,FunctionSpaceBase &
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();
@@ -86,6 +100,12 @@ template<class Assembler> void FixDofs(Assembler &assembler,std::vector<Dof> &do
}
class FilterDof
+{
+ public:
+ virtual bool operator()(Dof key)=0;
+};
+
+class FilterDofTrivial
{
public:
virtual bool operator()(Dof key) {return true;}
@@ -100,43 +120,60 @@ class FilterDofComponent :public FilterDof
{
int type=key.getType();
int icomp,iphys;
- Dof::getTwoIntsFromType(type, iphys, icomp);
+ Dof::getTwoIntsFromType(type, icomp, iphys);
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)
+template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MElement *e,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter)
{
- for (Iterator it=itbegin;it!=itend;++it)
+ std::vector<MVertex*> tabV;
+ int nv=e->getNumVertices();
+ std::vector<Dof> R;
+ 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;
- std::vector<MVertex*> tabV;
- int nv=e->getNumVertices();
- std::vector<Dof> R;
- 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)
+ Dof key=*itd;
+ if (filter(key))
{
- Dof key=*itd;
- if (filter(key))
+ for (int i=0;i<nv;++i)
{
- for (int i=0;i<nv;++i)
+ if (tabV[i]->getNum()==key.getEntity())
{
- if (tabV[i]->getNum()==key.getEntity())
- {
- assembler.fixDof(key, fct(tabV[i]->x(),tabV[i]->y(),tabV[i]->z()));
- break;
- }
+ assembler.fixDof(key, fct(tabV[i]->x(),tabV[i]->y(),tabV[i]->z()));
+ break;
}
}
}
}
}
+template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MVertex *v,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter)
+{
+ std::vector<Dof> R;
+ space.getKeys(v,R);
+ for (std::vector<Dof>::iterator itd=R.begin();itd!=R.end();++itd)
+ {
+ Dof key=*itd;
+ if (filter(key))
+ {
+ assembler.fixDof(key, fct(v->x(),v->y(),v->z()));
+ }
+ }
+}
+
+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)
+ {
+ FixNodalDofs(space,*it,assembler,fct,filter);
+ }
+}
template<class Iterator,class Assembler> void NumberDofs(FunctionSpaceBase &space,Iterator itbegin,Iterator itend,Assembler &assembler)
{
diff --git a/Solver/terms.h b/Solver/terms.h
index 721574d943..b5734be749 100644
--- a/Solver/terms.h
+++ b/Solver/terms.h
@@ -22,6 +22,7 @@
#include "groupOfElements.h"
+
// evaluation of a field ???
template<class T>
class Field {
@@ -80,6 +81,7 @@ class LinearTermBase
{
public:
virtual void get(MElement *ele,int npts,IntPt *GP,fullVector<double> &v) =0;
+ virtual void get(MVertex *v,fullVector<double> &m) =0;
};
template<class S1> class LinearTerm : public LinearTermBase
@@ -94,7 +96,7 @@ template<class S1> class LinearTerm : public LinearTermBase
class ScalarTermBase
{
public :
- virtual void get(MElement *ele,int npts,IntPt *GP,double &v) =0;
+ virtual void get(MElement *ele,int npts,IntPt *GP,double &val) =0;
};
class ScalarTerm : public ScalarTermBase
@@ -105,6 +107,60 @@ class ScalarTerm : public ScalarTermBase
+template<class S1,class S2> class LaplaceTerm : public BilinearTerm<S1,S2>
+{
+ public :
+ LaplaceTerm(S1& space1_,S2& space2_) : BilinearTerm<S1,S2>(space1_,space2_)
+ {}
+
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
+ {
+ Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
+ }
+ virtual void get(MVertex *v,fullMatrix<double> &m)
+ {
+ Msg::Error("LaplaceTerm<S1,S2> w/ S1!=S2 not implemented");
+ }
+}; // class
+
+
+
+template<class S1> class LaplaceTerm<S1,S1> : public BilinearTerm<S1,S1> // symmetric
+{
+ public :
+ LaplaceTerm(S1& space1_) : BilinearTerm<S1,S1>(space1_,space1_)
+ {}
+
+ virtual void get(MElement *ele,int npts,IntPt *GP,fullMatrix<double> &m)
+ {
+ int nbFF = BilinearTerm<S1,S1>::space1.getNumKeys(ele);
+ double jac[3][3];
+ m.resize(nbFF, nbFF);
+ m.setAll(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::GradType> Grads;
+ BilinearTerm<S1,S1>::space1.gradf(ele,u, v, w, Grads);
+ for (int j = 0; j < nbFF; j++)
+ {
+ for (int k = j; k < nbFF; k++)
+ {
+ double contrib=weight * detJ * dot(Grads[j],Grads[k]);
+ m(j,k)+=contrib;
+ if (j!=k) m(k,j)+=contrib;
+ }
+ }
+ }
+// m.print("");
+// exit(0);
+ }
+}; // class
+
+
+
+
template<class S1,class S2> class IsotropicElasticTerm : public BilinearTerm<S1,S2> // non symmetric
{
protected :
@@ -238,9 +294,9 @@ inline double dot(const double &a, const double &b)
template<class S1> class LoadTerm : public LinearTerm<S1>
{
- simpleFunction<typename S1::ValType> Load;
+ simpleFunction<typename S1::ValType> &Load;
public :
- LoadTerm(S1& space1_,simpleFunction<typename S1::ValType> Load_) :LinearTerm<S1>(space1_),Load(Load_) {};
+ LoadTerm(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);
@@ -253,12 +309,29 @@ template<class S1> class LoadTerm : public LinearTerm<S1>
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);
+ SPoint3 p;
+ ele->pnt(u, v, w, p);
+ typename S1::ValType load=Load(p.x(),p.y(),p.z());
for (int j = 0; j < nbFF ; ++j)
{
- m(j)+=dot(Vals[j],Load(u,v,w))*weight*detJ;
+ m(j)+=dot(Vals[j],load)*weight*detJ;
}
}
}
+
+ virtual void get(MVertex *ver,fullVector<double> &m)
+ {
+ double nbFF=LinearTerm<S1>::space1.getNumKeys(ver);
+ double jac[3][3];
+ m.resize(nbFF);
+ std::vector<typename S1::ValType> Vals;
+ LinearTerm<S1>::space1.f(ver, Vals);
+ typename S1::ValType load=Load(ver->x(),ver->y(),ver->z());
+ for (int j = 0; j < nbFF ; ++j)
+ {
+ m(j)=dot(Vals[j],load);
+ }
+ }
};
--
GitLab