Skip to content
Snippets Groups Projects
Commit 286a025f authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

*** empty log message ***

parent 6c1e0708
No related branches found
No related tags found
No related merge requests found
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#include <map> #include <map>
#include <vector> #include <vector>
#include "GmshMatrix.h" #include "GmshMatrix.h"
#include "gmshFunction.h"
#include "gmshAssembler.h" #include "gmshAssembler.h"
#include "GModel.h" #include "GModel.h"
#include "MElement.h" #include "MElement.h"
...@@ -97,63 +98,66 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> { ...@@ -97,63 +98,66 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
} }
} }
} }
void addDirichlet(int physical, void addDirichlet(int physical,
int dim, int dim,
int comp, int comp,
int field, int field,
const gmshFunction<scalar> & e, const gmshFunction<scalar> &e,
gmshAssembler<scalar> &lsys) gmshAssembler<scalar> &lsys)
{ {
std::vector<MVertex *> v; std::vector<MVertex *> v;
_gm->getMeshVertices (physical, dim,v); GModel *m = gmshTermOfFormulation<scalar>::_gm;
for (int i=0;i<v.size();i++) m->getMeshVertices(physical, dim, v);
lsys.fixVertex (v[i], comp, field , e(v[i]->x(),v[i]->y(),v[i]->z())); for (unsigned int i = 0; i < v.size(); i++)
lsys.fixVertex(v[i], comp, field, e(v[i]->x(), v[i]->y(), v[i]->z()));
} }
void addNeumann(int physical, void addNeumann(int physical,
int dim, int dim,
int comp, int comp,
int field, int field,
const gmshFunction<scalar> & fct, const gmshFunction<scalar> &fct,
gmshAssembler<scalar> &lsys) gmshAssembler<scalar> &lsys)
{ {
std::map<int, std::vector<GEntity*> > groups[4]; std::map<int, std::vector<GEntity*> > groups[4];
_gm->getPhysicalGroups(groups); GModel *m = gmshTermOfFormulation<scalar>::_gm;
std::map<int, std::vector<GEntity*> >::iterator it = groups[dim].find(physical); m->getPhysicalGroups(groups);
if (it == groups[dim].end())return; std::map<int, std::vector<GEntity*> >::iterator it = groups[dim].find(physical);
double jac[3][3]; if (it == groups[dim].end()) return;
double jac[3][3];
for (int i=0;i<it->second.size();++i){ for (unsigned int i = 0; i < it->second.size(); ++i){
GEntity *ge = it->second[i]; GEntity *ge = it->second[i];
for (int j=0; j<ge->getNumMeshElements();j++){ for (unsigned int j = 0; j < ge->getNumMeshElements(); j++){
MElement *e = ge->getMeshElement(j); MElement *e = ge->getMeshElement(j);
int integrationOrder = 2*e->getPolynomialOrder(); int integrationOrder = 2 * e->getPolynomialOrder();
int nbNodes = e->getNumVertices(); int nbNodes = e->getNumVertices();
int npts; int npts;
IntPt *GP; IntPt *GP;
e->getIntegrationPoints(integrationOrder, &npts, &GP); e->getIntegrationPoints(integrationOrder, &npts, &GP);
for (int ip=0;ip<npts;ip++){ for (int ip = 0; ip < npts; ip++){
const double u = GP[ip].pt[0]; const double u = GP[ip].pt[0];
const double v = GP[ip].pt[1]; const double v = GP[ip].pt[1];
const double w = GP[ip].pt[2]; const double w = GP[ip].pt[2];
const double weight = GP[ip].weight; const double weight = GP[ip].weight;
const double detJ = e->getJacobian (u,v,w,jac); const double detJ = e->getJacobian(u,v,w,jac);
double x = 0; double x = 0;
double y = 0; double y = 0;
double z = 0; double z = 0;
double sf[nbNodes]; double sf[nbNodes];
e->getShapeFunctions (u,v,w,sf); e->getShapeFunctions(u,v,w,sf);
for (int k=0;k<nbNodes;k++){ for (int k = 0; k < nbNodes; k++){
x += e->getVertex(k)->x() * sf[k]; x += e->getVertex(k)->x() * sf[k];
y += e->getVertex(k)->y() * sf[k]; y += e->getVertex(k)->y() * sf[k];
z += e->getVertex(k)->z() * sf[k]; z += e->getVertex(k)->z() * sf[k];
} }
const scalar FCT = fct (x,y,z); const scalar FCT = fct (x,y,z);
for (int k=0;k<nbNodes;k++){ for (int k = 0; k < nbNodes; k++){
lsys.assemble(e->getVertex(k), comp, field ,detJ * weight * sf[k] * FCT); lsys.assemble(e->getVertex(k), comp, field, detJ * weight * sf[k] * FCT);
} }
}
} }
} }
} }
}; };
#endif #endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment