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

* abstract elementary ("per element") solver data (mesh + the rest) in SElement 

* rewrite femTerm in terms of SElement
parent 5e5ec7ac
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Mesh/highOrderSmoother.cpp
+ 27
− 25
View file @ 0d376ca7
......@@ -115,12 +115,14 @@ struct p2data{
3 * t1->getNumVertices());
m2 = new fullMatrix<double>(3 * t2->getNumVertices(),
3 * t2->getNumVertices());
el.elementMatrix(t1,*m1);
el.elementMatrix(t2,*m2);
SElement se1(t1), se2(t2);
el.elementMatrix(&se1, *m1);
el.elementMatrix(&se2, *m2);
s->moveToTargetLocation(t1);
s->moveToTargetLocation(t2);
}
~p2data(){
~p2data()
{
delete m1;
delete m2;
}
......@@ -143,12 +145,14 @@ struct pNdata{
3 * t1->getNumVertices());
m2 = new fullMatrix<double>(3 * t2->getNumVertices(),
3 * t2->getNumVertices());
el.elementMatrix(t1,*m1);
el.elementMatrix(t2,*m2);
SElement se1(t1), se2(t2);
el.elementMatrix(&se1, *m1);
el.elementMatrix(&se2, *m2);
s->moveToTargetLocation(t1);
s->moveToTargetLocation(t2);
}
~pNdata(){
~pNdata()
{
delete m1;
delete m2;
}
......@@ -564,20 +568,17 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
fullVector<double> R2(n2);
fullMatrix<double> J23K33(n2, n3);
K33.set_all(0.0);
El.elementMatrix(e, K33);
SElement se(e);
El.elementMatrix(&se, K33);
computeMetricVector(gf, e, J32, J23, D3);
J23K33.gemm(J23, K33, 1, 0);
K22.gemm(J23K33, J32, 1, 0);
J23K33.mult(D3, R2);
for (int j = 0; j < n2; j++){
MVertex *vR;
int iCompR, iFieldR;
Dof RDOF = El.getLocalDofR(e, j);
Dof RDOF = El.getLocalDofR(&se, j);
myAssembler.assemble(RDOF, -R2(j));
for (int k = 0; k < n2; k++){
MVertex *vC;
int iCompC, iFieldC;
Dof CDOF = El.getLocalDofC(e, k);
Dof CDOF = El.getLocalDofC(&se, k);
myAssembler.assemble(RDOF, CDOF, K22(j, k));
}
}
......@@ -697,9 +698,10 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
// assembly of the elasticity term on the
// set of elements
for (int i=0;i<v.size();i++)
El.addToMatrix(myAssembler, v[i]);
for (int i = 0; i < v.size(); i++){
SElement se(v[i]);
El.addToMatrix(myAssembler, &se);
}
// solve the system
lsys->systemSolve();
}
......
Solver/SElement.h 0 → 100644
+ 26
− 0
View file @ 0d376ca7
// 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 _SELEMENT_H_
#define _SELEMENT_H_
#include "MElement.h"
// A solver element.
class SElement
{
private:
// the underlying mesh element
MElement *_e;
// store discrete function space and other data here
// ...
public:
SElement(MElement *e) : _e(e) {}
~SElement(){}
MElement *getMeshElement() const { return _e; }
};
#endif
Solver/elasticityTerm.cpp
+ 4
− 2
View file @ 0d376ca7
......@@ -6,8 +6,9 @@
#include "elasticityTerm.h"
#include "Numeric.h"
void elasticityTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
void elasticityTerm::elementMatrix(SElement *se, fullMatrix<double> &m) const
{
MElement *e = se->getMeshElement();
int nbNodes = e->getNumVertices();
int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
int npts;
......@@ -75,8 +76,9 @@ void elasticityTerm::elementMatrix(MElement *e, fullMatrix<double> &m) const
}
}
void elasticityTerm::elementVector(MElement *e, fullVector<double> &m) const
void elasticityTerm::elementVector(SElement *se, fullVector<double> &m) const
{
MElement *e = se->getMeshElement();
int nbNodes = e->getNumVertices();
int integrationOrder = 2 * e->getPolynomialOrder();
int npts;
......
Solver/elasticityTerm.h
+ 13
− 6
View file @ 0d376ca7
......@@ -9,7 +9,7 @@
#include "femTerm.h"
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "SElement.h"
#include "fullMatrix.h"
class elasticityTerm : public femTerm<double, double> {
......@@ -19,12 +19,19 @@ class elasticityTerm : public femTerm<double, double> {
SVector3 _volumeForce;
public:
// element matrix size : 3 dofs per vertex
virtual int sizeOfR(MElement *e) const { return 3 * e->getNumVertices(); }
virtual int sizeOfC(MElement *e) const { return 3 * e->getNumVertices(); }
virtual int sizeOfR(SElement *se) const
{
return 3 * se->getMeshElement()->getNumVertices();
}
virtual int sizeOfC(SElement *se) const
{
return 3 * se->getMeshElement()->getNumVertices();
}
// order dofs in the local matrix :
// dx1, dx2 ... dxn, dy1, ..., dyn, dz1, ... dzn
Dof getLocalDofR(MElement *e, int iRow) const
Dof getLocalDofR(SElement *se, int iRow) const
{
MElement *e = se->getMeshElement();
int iCompR = iRow / e->getNumVertices();
int ithLocalVertex = iRow % e->getNumVertices();
return Dof(e->getVertex(ithLocalVertex)->getNum(),
......@@ -34,8 +41,8 @@ class elasticityTerm : public femTerm<double, double> {
elasticityTerm(GModel *gm, double E, double nu, int iField) :
femTerm<double, double>(gm), _E(E), _nu(nu), _iField(iField) {}
void setVector(const SVector3 &f) { _volumeForce = f; }
void elementMatrix(MElement *e, fullMatrix<double> &m) const;
void elementVector(MElement *e, fullVector<double> &m) const;
void elementMatrix(SElement *se, fullMatrix<double> &m) const;
void elementVector(SElement *se, fullVector<double> &m) const;
};
#endif
Solver/femTerm.h
+ 23
− 23
View file @ 0d376ca7
......@@ -13,7 +13,7 @@
#include "simpleFunction.h"
#include "dofManager.h"
#include "GModel.h"
#include "MElement.h"
#include "SElement.h"
// a nodal finite element term : variables are always defined at nodes
// of the mesh
......@@ -25,20 +25,20 @@ class femTerm {
femTerm(GModel *gm) : _gm(gm) {}
virtual ~femTerm(){}
// return the number of columns of the element matrix
virtual int sizeOfC(MElement*) const = 0;
virtual int sizeOfC(SElement *se) const = 0;
// return the number of rows of the element matrix
virtual int sizeOfR(MElement*) const = 0;
virtual int sizeOfR(SElement *se) const = 0;
// in a given element, return the dof associated to a given row (column)
// of the local element matrix
virtual Dof getLocalDofR(MElement *e, int iRow) const = 0;
virtual Dof getLocalDofR(SElement *se, int iRow) const = 0;
// default behavior: symmetric
virtual Dof getLocalDofC(MElement *e, int iCol) const
virtual Dof getLocalDofC(SElement *se, int iCol) const
{
return getLocalDofR(e, iCol);
return getLocalDofR(se, iCol);
}
// compute the elementary matrix
virtual void elementMatrix(MElement *e, fullMatrix<dataMat> &m) const = 0;
virtual void elementVector(MElement *e, fullVector<dataVec> &m) const
virtual void elementMatrix(SElement *se, fullMatrix<dataMat> &m) const = 0;
virtual void elementVector(SElement *se, fullVector<dataVec> &m) const
{
m.scale(0.0);
}
......@@ -47,29 +47,29 @@ class femTerm {
void addToMatrix(dofManager<dataVec, dataMat> &dm, GEntity *ge) const
{
for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
MElement *e = ge->getMeshElement(i);
addToMatrix(dm, e);
SElement se(ge->getMeshElement(i));
addToMatrix(dm, &se);
}
}
// add the contribution from a single element to the dof manager
void addToMatrix(dofManager<dataVec, dataMat> &dm, MElement *e) const
void addToMatrix(dofManager<dataVec, dataMat> &dm, SElement *se) const
{
const int nbR = sizeOfR(e);
const int nbC = sizeOfC(e);
const int nbR = sizeOfR(se);
const int nbC = sizeOfC(se);
fullMatrix<dataMat> localMatrix(nbR, nbC);
elementMatrix(e, localMatrix);
addToMatrix(dm, localMatrix, e);
elementMatrix(se, localMatrix);
addToMatrix(dm, localMatrix, se);
}
void addToMatrix(dofManager<dataVec, dataMat> &dm,
fullMatrix<dataMat> &localMatrix,
MElement *e) const
SElement *se) const
{
const int nbR = localMatrix.size1();
const int nbC = localMatrix.size2();
for (int j = 0; j < nbR; j++){
Dof R = getLocalDofR(e, j);
Dof R = getLocalDofR(se, j);
for (int k = 0; k < nbC; k++){
Dof C = getLocalDofC(e, k);
Dof C = getLocalDofC(se, k);
dm.assemble(R, C, localMatrix(j, k));
}
}
......@@ -123,12 +123,12 @@ class femTerm {
void addToRightHandSide(dofManager<dataVec, dataMat> &dm, GEntity *ge) const
{
for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
MElement *e = ge->getMeshElement(i);
int nbR = sizeOfR(e);
SElement se(ge->getMeshElement(i));
int nbR = sizeOfR(&se);
fullVector<dataVec> V(nbR);
elementVector(e, V);
elementVector(&se, V);
// assembly
for (int j = 0; j < nbR; j++) dm.assemble(getLocalDofR(e, j), V(j));
for (int j = 0; j < nbR; j++) dm.assemble(getLocalDofR(&se, j), V(j));
}
}
};
......
Solver/helmholtzTerm.h
+ 15
− 8
View file @ 0d376ca7
......@@ -11,7 +11,7 @@
#include "simpleFunction.h"
#include "Gmsh.h"
#include "GModel.h"
#include "MElement.h"
#include "SElement.h"
#include "fullMatrix.h"
#include "Numeric.h"
......@@ -28,18 +28,25 @@ class helmoltzTerm : public femTerm<scalar, scalar> {
: femTerm<scalar,scalar>(gm), _k(k), _a(a), _iFieldR(iFieldR),
_iFieldC(iFieldC) {}
// one dof per vertex (nodal fem)
virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
Dof getLocalDofR(MElement *e, int iRow) const
virtual int sizeOfR(SElement *se) const
{
return Dof(e->getVertex(iRow)->getNum(), _iFieldR);
return se->getMeshElement()->getNumVertices();
}
Dof getLocalDofC(MElement *e, int iRow) const
virtual int sizeOfC(SElement *se) const
{
return Dof(e->getVertex(iRow)->getNum(), _iFieldC);
return se->getMeshElement()->getNumVertices();
}
virtual void elementMatrix(MElement *e, fullMatrix<scalar> &m) const
Dof getLocalDofR(SElement *se, int iRow) const
{
return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldR);
}
Dof getLocalDofC(SElement *se, int iRow) const
{
return Dof(se->getMeshElement()->getVertex(iRow)->getNum(), _iFieldC);
}
virtual void elementMatrix(SElement *se, fullMatrix<scalar> &m) const
{
MElement *e = se->getMeshElement();
// compute integration rule
const int integrationOrder = (_a) ? 2 * e->getPolynomialOrder() :
2 * (e->getPolynomialOrder() - 1);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment