Skip to content
Snippets Groups Projects
Commit 6146be39 authored by Éric Béchet's avatar Éric Béchet
Browse files

weight watchers

parent fd1b5772
No related branches found
No related tags found
No related merge requests found
......@@ -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);
......
......@@ -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);
}
}
/*
printf("-- done assembling!\n");
for (int i=0;i<pAssembler->sizeOfR();i++)
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); i++)
{
if (fabs(lsys->getFromRightHandSide(i))>0.000001)
printf("%g ",lsys->getFromRightHandSide(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("\n");
*/
// solving
printf("-- done assembling!\n");
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);
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);
std::map<int, std::vector<GEntity*> > groups[4];
pModel->getPhysicalGroups(groups);
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];
// 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);
double s[9] = {sxx, sxy, sxz,
sxy, syy, syz,
sxz, syz, szz};
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);
return ComputeVonMises(s);
}
for (unsigned int i = 0; i < edgeDiri.size(); i++)
PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
{
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 (unsigned int i = 0; i < faceDiri.size(); i++)
std::set<MVertex*> v;
for (unsigned int i = 0; i < elasticFields.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);
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));
}
// we number the nodes : when a node is numbered, it cannot be numbered
// again with another number. so we loop over elements
}
std::map<int, std::vector<double> > data;
for ( std::set<MVertex*>::iterator it = v.begin(); it != v.end(); ++it)
{
SVector3 val(0,0,0);
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());
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;
}
GaussQuadrature Integ_Boundary(GaussQuadrature::Val);
for (unsigned int i = 0; i < edgeNeu.size(); i++)
PView *elasticitySolver::buildVonMisesView(const std::string &postFileName)
{
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);
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 < 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);
}
for (unsigned int i = 0; i < volumeNeu.size(); i++)
#else
PView* elasticitySolver::buildDisplacementView (const std::string &postFileName)
{
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);
Msg::Error("Post-pro module not available");
return 0;
}
GaussQuadrature Integ_Bulk(GaussQuadrature::GradGrad);
for (unsigned int i = 0; i < elasticFields.size(); i++)
PView* elasticitySolver::buildVonMisesView(const std::string &postFileName)
{
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);
Msg::Error("Post-pro module not available");
return 0;
}
printf("-- done assembling!\n");
lsys->systemSolve();
printf("-- done solving!\n");
}
#endif
......@@ -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
......@@ -68,12 +68,14 @@ 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 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
......@@ -90,25 +93,23 @@ class FunctionSpace : public FunctionSpaceBase
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 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;};
......@@ -119,6 +120,7 @@ class ScalarLagrangeFunctionSpace : public FunctionSpace<double>
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);
};
......
......@@ -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,18 +120,15 @@ 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)
{
for (Iterator it=itbegin;it!=itend;++it)
template<class Assembler> void FixNodalDofs(FunctionSpaceBase &space,MElement *e,Assembler &assembler,simpleFunction<typename Assembler::dataVec> &fct,FilterDof &filter)
{
MElement *e=*it;
std::vector<MVertex*> tabV;
int nv=e->getNumVertices();
std::vector<Dof> R;
......@@ -135,8 +152,28 @@ template<class Iterator,class Assembler> void FixNodalDofs(FunctionSpaceBase &sp
}
}
}
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)
{
......
......@@ -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);
}
}
};
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment