From 6dd5c4d3ac468aac16716e6b907124baafaaa32e Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Mon, 23 Mar 2009 15:30:01 +0000
Subject: [PATCH] *** empty log message ***
---
Makefile | 2 +-
Numeric/gmshElasticity.cpp | 26 ++++++++++++++++++++++++++
Numeric/gmshElasticity.h | 7 +++++--
Numeric/gmshTermOfFormulation.h | 31 +++++++++++++++++++++++++++++++
Post/PViewDataGModelIO.cpp | 1 -
5 files changed, 63 insertions(+), 4 deletions(-)
diff --git a/Makefile b/Makefile
index 9e24091264..46059a80a8 100644
--- a/Makefile
+++ b/Makefile
@@ -20,7 +20,7 @@ GMSH_DATE = `date "+%Y%m%d"`
GMSH_API = Geo/GModel.h Geo/GEntity.h Geo/GPoint.h\
Geo/GVertex.h Geo/GEdge.h Geo/GFace.h Geo/GRegion.h\
Geo/GEdgeLoop.h Geo/GFaceCompound.h\
- Geo/MVertex.h Geo/MEdge.h Geo/MFace.h Geo/MElement.h\
+ Geo/MVertex.h Geo/MEdge.h Geo/MFace.h Geo/MElement.h Geo/MTriangle.h Geo/MQuadrangle.h Geo/MTetrahedron.h Geo/MHexahedron.h Geo/MPyramid.h Geo/MPrism.h Geo/MLine.h \
Geo/discreteVertex.h Geo/discreteEdge.h Geo/discreteFace.h Geo/discreteRegion.h\
Geo/SPoint2.h Geo/SPoint3.h Geo/SVector3.h Geo/SBoundingBox3d.h\
Geo/Pair.h Geo/Range.h\
diff --git a/Numeric/gmshElasticity.cpp b/Numeric/gmshElasticity.cpp
index 379e72d64b..2750a47ff1 100644
--- a/Numeric/gmshElasticity.cpp
+++ b/Numeric/gmshElasticity.cpp
@@ -74,3 +74,29 @@ void gmshElasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
m.gemm(BTH, B, weight * detJ, 1.);
}
}
+
+void gmshElasticityTerm::elementVector(MElement *e, gmshVector<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;
+ }
+ }
+}
diff --git a/Numeric/gmshElasticity.h b/Numeric/gmshElasticity.h
index 348545483a..037aeedbd2 100644
--- a/Numeric/gmshElasticity.h
+++ b/Numeric/gmshElasticity.h
@@ -15,6 +15,7 @@
class gmshElasticityTerm : public gmshNodalFemTerm<double> {
double _E, _nu;
int _iField;
+ SVector3 _f;
protected:
virtual int sizeOfR(MElement *e) const { return 3 * e->getNumVertices(); }
virtual int sizeOfC(MElement *e) const { return 3 * e->getNumVertices(); }
@@ -26,9 +27,11 @@ class gmshElasticityTerm : public gmshNodalFemTerm<double> {
*iFieldR = _iField;
}
public:
- gmshElasticityTerm(GModel *gm, double E, double nu, int iField = 1) :
- gmshNodalFemTerm<double>(gm), _E(E), _nu(nu), _iField(iField){}
+ 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, gmshMatrix<double> &m) const;
+ void elementVector(MElement *e, gmshVector<double> &m) const;
};
#endif
diff --git a/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h
index fd264ef7bc..2ee14a8d3c 100644
--- a/Numeric/gmshTermOfFormulation.h
+++ b/Numeric/gmshTermOfFormulation.h
@@ -47,6 +47,9 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
gmshNodalFemTerm(GModel *gm) : gmshTermOfFormulation<scalar>(gm) {}
virtual ~gmshNodalFemTerm (){}
virtual void elementMatrix(MElement *e, gmshMatrix<scalar> &m) const = 0;
+ virtual void elementVector(MElement *e, gmshVector<scalar> &m) const {
+ m.scale(0.0);
+ }
void addToMatrix(gmshAssembler<scalar> &lsys) const
{
GModel *m = gmshTermOfFormulation<scalar>::_gm;
@@ -151,6 +154,34 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
}
}
}
+ void addToRightHandSide (gmshAssembler<scalar> &lsys) const {
+ GModel *m = gmshTermOfFormulation<scalar>::_gm;
+ if (m->getNumRegions()){
+ for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
+ addToRightHandSide(lsys,*it);
+ }
+ }
+ else if(m->getNumFaces()){
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+ addToRightHandSide(lsys,*it);
+ }
+ }
+ }
+ void addToRightHandSide (gmshAssembler<scalar> &lsys, GEntity *ge) const {
+ for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
+ MElement *e = ge->getMeshElement (i);
+ int nbR = sizeOfR(e);
+ gmshVector<scalar> V (nbR);
+ elementVector (e, V);
+ // assembly
+ for (int j=0;j<nbR;j++){
+ MVertex *vR;int iCompR,iFieldR;
+ getLocalDofR (e,j,&vR,&iCompR,&iFieldR);
+ lsys.assemble(vR,iCompR,iFieldR,V(j));
+ }
+ }
+ }
+
};
#endif
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 9a0344dd5b..9ad02b24bb 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -120,7 +120,6 @@ bool PViewDataGModel::writeMSH(std::string fileName, bool binary)
GModel *model = _steps[0]->getModel();
binary = true;
-
if(!model->writeMSH(fileName, 2.0, binary, true)) return false;
// append data
--
GitLab