Skip to content
Snippets Groups Projects
Commit c0d20d98 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

*** empty log message ***

parent 04a6e451
No related branches found
No related tags found
No related merge requests found
......@@ -70,9 +70,9 @@ set(GMSH_API
Mesh/meshGFaceDelaunayInsertion.h
Post/PView.h Post/PViewData.h Plugin/PluginManager.h
Graphics/drawContext.h
Solver/dofManager.h Solver/elasticityTerm.h Solver/femTerm.h
Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
Numeric/gmshAssembler.h Numeric/gmshTermOfFormulation.h
Numeric/gmshLaplace.h Numeric/gmshElasticity.h Numeric/gmshLinearSystem.h
Numeric/gmshLaplace.h Numeric/gmshLinearSystem.h
Numeric/gmshLinearSystemGmm.h Numeric/gmshLinearSystemFull.h
Numeric/gmshLinearSystemCSR.h
contrib/kbipack/gmp_normal_form.h contrib/kbipack/gmp_matrix.h
......
......@@ -21,7 +21,6 @@ set(SRC
gmshCrossConf.cpp
gmshDistance.cpp
gmshConvexCombination.cpp
gmshElasticity.cpp
gmshLinearSystemCSR.cpp
)
......
// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include "gmshElasticity.h"
#include "Numeric.h"
/*
Non linear version :
-) \sigma_{}
*/
void gmshElasticityTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
{
int nbNodes = e->getNumVertices();
int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
int npts;
IntPt *GP;
double jac[3][3];
double invjac[3][3];
double Grads[256][3], grads[100][3];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
m.set_all(0.);
double FACT = _E / (1 + _nu);
double C11 = FACT * (1 - _nu) / (1 - 2 * _nu);
double C12 = FACT * _nu / (1 - 2 * _nu);
double C44 = (C11 - C12) / 2;
const double C[6][6] =
{ {C11, C12, C12, 0, 0, 0},
{C12, C11, C12, 0, 0, 0},
{C12, C12, C11, 0, 0, 0},
{ 0, 0, 0, C44, 0, 0},
{ 0, 0, 0, 0, C44, 0},
{ 0, 0, 0, 0, 0, C44} };
fullMatrix<double> H(6, 6);
fullMatrix<double> B(6, 3 * nbNodes);
fullMatrix<double> BTH(3 * nbNodes, 6);
fullMatrix<double> BT(3 * nbNodes, 6);
for (int i = 0; i < 6; i++)
for (int j = 0; j < 6; j++)
H(i, j) = C[i][j];
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 = e->getJacobian(u, v, w, jac);
e->getGradShapeFunctions(u, v, w, grads);
inv3x3(jac, invjac);
B.set_all(0.);
BT.set_all(0.);
for (int j = 0; j < nbNodes; j++){
Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
invjac[0][2] * grads[j][2];
Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
invjac[1][2] * grads[j][2];
Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
invjac[2][2] * grads[j][2];
BT(j, 0) = B(0, j) = Grads[j][0];
BT(j, 3) = B(3, j) = Grads[j][1];
BT(j, 4) = B(4, j) = Grads[j][2];
BT(j + nbNodes, 1) = B(1, j + nbNodes) = Grads[j][1];
BT(j + nbNodes, 3) = B(3, j + nbNodes) = Grads[j][0];
BT(j + nbNodes, 5) = B(5, j + nbNodes) = Grads[j][2];
BT(j + 2 * nbNodes, 2) = B(2, j + 2 * nbNodes) = Grads[j][2];
BT(j + 2 * nbNodes, 4) = B(4, j + 2 * nbNodes) = Grads[j][0];
BT(j + 2 * nbNodes, 5) = B(5, j + 2 * nbNodes) = Grads[j][1];
}
BTH.set_all(0.);
BTH.gemm(BT, H);
m.gemm(BTH, B, weight * detJ, 1.);
}
}
void gmshElasticityTerm::elementVector(MElement *e, fullVector<double> &m) const{
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * e->getPolynomialOrder();
int npts;
IntPt *GP;
double jac[3][3];
double ff[256];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
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 = e->getJacobian(u, v, w, jac);
e->getShapeFunctions(u, v, w, ff);
for (int j = 0; j < nbNodes; j++){
m(j) += _f.x() *ff[j] * weight * detJ *.5;
m(j+nbNodes) += _f.y() *ff[j] * weight * detJ *.5;
m(j+2*nbNodes) += _f.z() *ff[j] * weight * detJ *.5;
}
}
}
// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#ifndef _GMSH_ELASTICITY_H_
#define _GMSH_ELASTICITY_H_
#include "gmshTermOfFormulation.h"
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "fullMatrix.h"
class gmshElasticityTerm : public gmshNodalFemTerm<double> {
protected:
double _E, _nu;
int _iField;
SVector3 _f;
public:
virtual int sizeOfR(MElement *e) const { return 3 * e->getNumVertices(); }
virtual int sizeOfC(MElement *e) const { return 3 * e->getNumVertices(); }
void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
{
*iCompR = iRow / e->getNumVertices();
int ithLocalVertex = iRow % e->getNumVertices();
*vR = e->getVertex(ithLocalVertex);
*iFieldR = _iField;
}
public:
gmshElasticityTerm(GModel *gm, double E, double nu, int iField = 1) :
gmshNodalFemTerm<double>(gm), _E(E), _nu(nu), _iField(iField){}
void setVector(const SVector3 &f) {_f = f;}
void elementMatrix(MElement *e, fullMatrix<double> &m) const;
void elementVector(MElement *e, fullVector<double> &m) const;
};
#endif
// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include "gmshHelmholtz.h"
#include "Numeric.h"
void gmshHelmholtzTerm::elementMatrix(MElement *e,
fullMatrix<std::complex<double> > &m) const
{
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * e->getPolynomialOrder();
int npts;
IntPt *GP;
double jac[3][3];
double invjac[3][3];
double sf[256], Grads[256][3], grads[256][3];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
m.set_all(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 = e->getJacobian(u, v, w, jac);
SPoint3 p; e->pnt(u, v, w, p);
const std::complex<double> kk = (*_waveNumber)(p.x(), p.y(), p.z());
inv3x3(jac, invjac) ;
e->getShapeFunctions(u, v, w, sf);
e->getGradShapeFunctions(u, v, w, grads);
for (int j = 0; j < nbNodes; j++){
Grads[j][0] = invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
invjac[0][2] * grads[j][2];
Grads[j][1] = invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
invjac[1][2] * grads[j][2];
Grads[j][2] = invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
invjac[2][2] * grads[j][2];
}
for (int j = 0; j < nbNodes; j++){
for (int k = 0; k < nbNodes; k++){
m(j, k) += ((Grads[j][0] * Grads[k][0] +
Grads[j][1] * Grads[k][1] +
Grads[j][2] * Grads[k][2])
- kk * kk * sf[j] * sf[k]) * weight * detJ;
}
}
}
}
// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#ifndef _GMSH_HELMHOLTZ_H_
#define _GMSH_HELMHOLTZ_H_
#include <complex>
#include "gmshTermOfFormulation.h"
#include "simpleFunction.h"
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "fullMatrix.h"
class gmshHelmholtzTerm : public gmshNodalFemTerm<std::complex<double> > {
private:
const simpleFunction<std::complex<double> > *_waveNumber;
const int _iField ;
protected:
virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
{
*vR = e->getVertex(iRow);
*iCompR = 0;
*iFieldR = _iField;
}
public:
gmshHelmholtzTerm(GModel *gm, simpleFunction<std::complex<double> > *waveNumber,
int iField = 0)
: gmshNodalFemTerm<std::complex<double> >(gm), _waveNumber(waveNumber), _iField(iField){}
virtual void elementMatrix(MElement *e, fullMatrix<std::complex<double> > &m) const;
};
#endif
// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include "gmshProjection.h"
#include "Numeric.h"
void gmshProjectionTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
{
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * e->getPolynomialOrder();
int npts;
IntPt *GP;
double jac[3][3];
double sf[256];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
m.set_all(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 = e->getJacobian(u, v, w, jac);
SPoint3 p; e->pnt(u, v, w, p);
e->getShapeFunctions(u, v, w, sf);
for (int j = 0; j < nbNodes; j++){
for (int k = 0; k < nbNodes; k++){
m(j, k) += sf[j] * sf[k] * weight * detJ;
}
}
}
}
// THIS FILE WILL BE REMOVED AS SOON AS THE NEW SOLVER INTERFACE IS READY
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#ifndef _GMSH_PROJECTION_H_
#define _GMSH_PROJECTION_H_
#include "gmshTermOfFormulation.h"
#include "simpleFunction.h"
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "fullMatrix.h"
class gmshProjectionTerm : public gmshNodalFemTerm<double> {
private:
const int _iField ;
protected:
virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
{
*vR = e->getVertex(iRow);
*iCompR = 0;
*iFieldR = _iField;
}
public:
gmshProjectionTerm(GModel *gm, int iField = 0) :
gmshNodalFemTerm<double>(gm), _iField(iField){}
virtual void elementMatrix(MElement *e, fullMatrix<double> &m) const;
};
#endif
......@@ -23,7 +23,6 @@ set(SRC
Annotate.cpp Remove.cpp
Probe.cpp
HarmonicToTime.cpp ModulusPhase.cpp
FiniteElement.cpp
HomologyComputation.cpp
)
......
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#include <stdlib.h>
#include "GmshConfig.h"
#include "GModel.h"
#include "FiniteElement.h"
#include "gmshAssembler.h"
#include "gmshProjection.h"
#include "gmshLaplace.h"
#include "gmshHelmholtz.h"
#include "gmshLinearSystemGmm.h"
#include "mathEvaluator.h"
StringXNumber FiniteElementOptions_Number[] = {
{GMSH_FULLRC, "Omega", NULL, 1.},
{GMSH_FULLRC, "Gamma1", NULL, 0.},
{GMSH_FULLRC, "Gamma1Value", NULL, 0.},
{GMSH_FULLRC, "Gamma2", NULL, 0.},
{GMSH_FULLRC, "Gamma2Value", NULL, 0.}
};
StringXString FiniteElementOptions_String[] = {
{GMSH_FULLRC, "Equation", NULL, "Projection"},
{GMSH_FULLRC, "Parameter", NULL, "Sin(x*y)"},
{GMSH_FULLRC, "Gamma1BC", NULL, ""},
{GMSH_FULLRC, "Gamma2BC", NULL, ""},
};
extern "C"
{
GMSH_Plugin *GMSH_RegisterFiniteElementPlugin()
{
return new GMSH_FiniteElementPlugin();
}
}
std::string GMSH_FiniteElementPlugin::getHelp() const
{
return "Plugin(FiniteElement) solves simple PDEs\n"
"using the finite element method. This is only\n"
"intended as a demonstration tool: it is NOT\n"
"intended for general use."
"\n"
"Plugin(FiniteElement) creates a new view.\n";
}
int GMSH_FiniteElementPlugin::getNbOptions() const
{
return sizeof(FiniteElementOptions_Number) / sizeof(StringXNumber);
}
StringXNumber *GMSH_FiniteElementPlugin::getOption(int iopt)
{
return &FiniteElementOptions_Number[iopt];
}
int GMSH_FiniteElementPlugin::getNbOptionsStr() const
{
return sizeof(FiniteElementOptions_String) / sizeof(StringXString);
}
StringXString *GMSH_FiniteElementPlugin::getOptionStr(int iopt)
{
return &FiniteElementOptions_String[iopt];
}
template<class scalar>
class gmshMathEvalFunction : public simpleFunction<scalar> {
private:
std::string _str;
mathEvaluator *_f;
public:
gmshMathEvalFunction(std::string str) : _str(str)
{
std::vector<std::string> expressions(1), variables(3);
expressions[0] = str;
variables[0] = "x";
variables[1] = "y";
variables[2] = "z";
_f = new mathEvaluator(expressions, variables);
if(expressions.empty()){
delete _f;
_f = 0;
}
}
virtual ~gmshMathEvalFunction()
{
if(_f) delete _f;
}
virtual scalar operator () (double x, double y, double z) const
{
if(_f){
std::vector<double> values(3), res(1);
values[0] = x;
values[1] = y;
values[2] = z;
if(_f->eval(values, res)) return res[0];
}
return atof(_str.c_str());
}
};
#if defined(HAVE_GMM)
template<class scalar>
class solver{
public:
gmshLinearSystemGmm<scalar> *lsys;
gmshAssembler<scalar> *myAssembler;
solver()
{
int volume = (int)FiniteElementOptions_Number[0].def;
int surface1 = (int)FiniteElementOptions_Number[1].def;
double value1 = FiniteElementOptions_Number[2].def;
int surface2 = (int)FiniteElementOptions_Number[3].def;
double value2 = FiniteElementOptions_Number[4].def;
std::string bc1 = FiniteElementOptions_String[2].def;
std::string bc2 = FiniteElementOptions_String[3].def;
lsys = new gmshLinearSystemGmm<scalar>;
lsys->setPrec(1.e-10);
lsys->setNoisy(2);
myAssembler = new gmshAssembler<scalar>(lsys);
GModel *m = GModel::current();
std::vector<MVertex*> vertices;
int dim = m->getNumRegions() ? 3 : 2;
if(bc1 == "Dirichlet"){
m->getMeshVerticesForPhysicalGroup(dim - 1, surface1, vertices);
for(unsigned int i = 0; i < vertices.size(); i++)
myAssembler->fixVertex(vertices[i], 0, 1, value1);
}
if(bc2 == "Dirichlet"){
m->getMeshVerticesForPhysicalGroup(dim - 1, surface2, vertices);
for(unsigned int i = 0; i < vertices.size(); i++)
myAssembler->fixVertex(vertices[i], 0, 1, value2);
}
m->getMeshVerticesForPhysicalGroup(dim, volume, vertices);
for(unsigned int i = 0; i < vertices.size(); i++)
myAssembler->numberVertex(vertices[i], 0, 1);
Msg::Info("Assembler: %d dofs, %d fixed",
myAssembler->sizeOfR(), myAssembler->sizeOfF());
}
};
#endif
PView *GMSH_FiniteElementPlugin::execute(PView *v)
{
std::string equation = FiniteElementOptions_String[0].def;
std::string parameter = FiniteElementOptions_String[1].def;
int volume = (int)FiniteElementOptions_Number[0].def;
#if defined(HAVE_GMM)
GModel *m = GModel::current();
std::map<int, std::vector<GEntity*> > groups[4];
m->getPhysicalGroups(groups);
int dim = m->getNumRegions() ? 3 : 2;
std::vector<MVertex*> vertices;
m->getMeshVerticesForPhysicalGroup(dim, volume, vertices);
std::map<int, std::vector<double> > data;
if(equation == "Projection"){
solver<double> s;
if(!s.myAssembler->sizeOfR()) return 0;
gmshMathEvalFunction<double> par(parameter);
gmshProjectionTerm projection(m, 1);
for(unsigned int i = 0; i < groups[dim][volume].size(); i++)
projection.addToMatrix(*s.myAssembler, groups[dim][volume][i]);
projection.addNeumann(volume, dim, 0, 1, par, *s.myAssembler);
s.lsys->systemSolve();
for (unsigned int i = 0; i < vertices.size(); i++)
data[vertices[i]->getNum()].push_back
(s.myAssembler->getDofValue(vertices[i], 0, 1));
new PView("projection", "NodeData", m, data, 0.);
}
else if(equation == "Laplace"){
solver<double> s;
if(!s.myAssembler->sizeOfR()) return 0;
gmshMathEvalFunction<double> par(parameter);
gmshLaplaceTerm laplace(m, &par, 1);
for(unsigned int i = 0; i < groups[dim][volume].size(); i++)
laplace.addToMatrix(*s.myAssembler, groups[dim][volume][i]);
s.lsys->systemSolve();
for (unsigned int i = 0; i < vertices.size(); i++)
data[vertices[i]->getNum()].push_back
(s.myAssembler->getDofValue(vertices[i], 0, 1));
new PView("laplace", "NodeData", m, data, 0.);
}
else if(equation == "Helmholtz"){
solver<std::complex<double> > s;
if(!s.myAssembler->sizeOfR()) return 0;
gmshMathEvalFunction<std::complex<double> > par(parameter);
gmshHelmholtzTerm helmholtz(m, &par, 1);
for(unsigned int i = 0; i < groups[dim][volume].size(); i++)
helmholtz.addToMatrix(*s.myAssembler, groups[dim][volume][i]);
s.lsys->setGmres(1);
s.lsys->systemSolve();
for (unsigned int i = 0; i < vertices.size(); i++)
data[vertices[i]->getNum()].push_back
(s.myAssembler->getDofValue(vertices[i], 0, 1).real());
PView *pv = new PView("helmholtz", "NodeData", m, data, 0.);
data.clear();
for (unsigned int i = 0; i < vertices.size(); i++)
data[vertices[i]->getNum()].push_back
(s.myAssembler->getDofValue(vertices[i], 0, 1).imag());
pv->addStep(m, data, 1.);
}
else{
Msg::Error("Unknown equation to solve");
}
#endif
return 0;
}
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
#ifndef _FINITE_ELEMENT_H_
#define _FINITE_ELEMENT_H_
#include "Plugin.h"
extern "C"
{
GMSH_Plugin *GMSH_RegisterFiniteElementPlugin();
}
class GMSH_FiniteElementPlugin : public GMSH_PostPlugin
{
public:
GMSH_FiniteElementPlugin(){}
std::string getName() const { return "FiniteElement"; }
std::string getHelp() const;
int getNbOptions() const;
StringXNumber* getOption(int iopt);
int getNbOptionsStr() const;
StringXString* getOptionStr(int iopt);
PView *execute(PView *);
};
#endif
......@@ -43,7 +43,6 @@
#include "Probe.h"
#include "FieldView.h"
#include "GSHHS.h"
#include "FiniteElement.h"
#include "HomologyComputation.h"
// for testing purposes only :-)
......@@ -216,8 +215,6 @@ void PluginManager::registerDefaultPlugins()
("Evaluate", GMSH_RegisterEvaluatePlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("CutParametric", GMSH_RegisterCutParametricPlugin()));
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("FiniteElement", GMSH_RegisterFiniteElementPlugin()));
#if defined(HAVE_KBIPACK)
allPlugins.insert(std::pair<std::string, GMSH_Plugin*>
("HomologyComputation", GMSH_RegisterHomologyComputationPlugin()));
......
......@@ -57,9 +57,7 @@ class elasticitySolver{
// std::pair<PView *, PView*> buildErrorEstimateView
// (const std::string &errorFileName, double, int);
// std::pair<PView *, PView*> buildErrorEstimateView
// (const std::string &errorFileName, const gmshElasticityData &ref, double, int);
// (const std::string &errorFileName, const elasticityData &ref, double, int);
};
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment