From 9bf05e3921a7907a165a2a6d1c3f5b6e9df9f727 Mon Sep 17 00:00:00 2001
From: Samuel Melchior <samuel.melchior@uclouvain.be>
Date: Tue, 8 Dec 2009 07:55:31 +0000
Subject: [PATCH] 1D
---
CMakeLists.txt | 2 +
Geo/MElement.cpp | 9 +++--
Geo/MElement.h | 3 ++
Geo/MLine.h | 1 +
Geo/MPoint.h | 12 ++++++
Numeric/polynomialBasis.cpp | 24 ++++++++++--
Numeric/polynomialBasis.h | 5 +++
Solver/dgAlgorithm.cpp | 5 +--
Solver/dgGroupOfElements.cpp | 73 ++++++++++++++++++++++++++++++------
Solver/dgGroupOfElements.h | 2 +
10 files changed, 115 insertions(+), 21 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 716c8f82b8..8fca86e060 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1016,3 +1016,5 @@ target_link_libraries(gmshdgsw ${LINK_LIBRARIES})
add_executable(gmshdgwave EXCLUDE_FROM_ALL Solver/dgMainWaveEquation.cpp ${GMSH_SRC})
target_link_libraries(gmshdgwave ${LINK_LIBRARIES})
+add_executable(gmshdgsam EXCLUDE_FROM_ALL Solver/dgMainSam.cpp ${GMSH_SRC})
+target_link_libraries(gmshdgsam ${LINK_LIBRARIES})
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index ca3ba8a384..2d475aaff1 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -154,6 +154,7 @@ static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
case 0:
{
+ dJ = 1.0;
jac[0][0] = jac[1][1] = jac[2][2] = 1.0;
jac[0][1] = jac[1][0] = jac[2][0] = 0.0;
jac[0][2] = jac[1][2] = jac[2][1] = 0.0;
@@ -165,9 +166,9 @@ static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
// regularize matrix
double a[3], b[3], c[3];
- a[0] = ele->getVertex(1)->x() - ele->getVertex(0)->x();
- a[1] = ele->getVertex(1)->y() - ele->getVertex(0)->y();
- a[2] = ele->getVertex(1)->z() - ele->getVertex(0)->z();
+ a[0] = jac[0][0];
+ a[1] = jac[0][1];
+ a[2] = jac[0][2];
if((fabs(a[0]) >= fabs(a[1]) && fabs(a[0]) >= fabs(a[2])) ||
(fabs(a[1]) >= fabs(a[0]) && fabs(a[1]) >= fabs(a[2]))) {
b[0] = a[1]; b[1] = -a[0]; b[2] = 0.;
@@ -175,7 +176,9 @@ static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
else {
b[0] = 0.; b[1] = a[2]; b[2] = -a[1];
}
+ norme(b);
prodve(a, b, c);
+ norme(c);
jac[0][1] = b[0]; jac[1][1] = b[1]; jac[2][1] = b[2];
jac[0][2] = c[0]; jac[1][2] = c[1]; jac[2][2] = c[2];
break;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 44f6c2b73c..90de79149e 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -73,6 +73,9 @@ class MElement
// get the vertices
virtual int getNumVertices() const = 0;
virtual MVertex *getVertex(int num) = 0;
+ // give an MVertex as input and get its local number
+ virtual void getVertexInfo (const MVertex * vertex, int &ithVertex) const
+ {throw;}
// get the vertex using the I-deas UNV ordering
virtual MVertex *getVertexUNV(int num){ return getVertex(num); }
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 74542824aa..7266b32e38 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -37,6 +37,7 @@ class MLine : public MElement {
virtual int getDim(){ return 1; }
virtual int getNumVertices() const { return 2; }
virtual MVertex *getVertex(int num){ return _v[num]; }
+ virtual void getVertexInfo (const MVertex * vertex, int &ithVertex) const { ithVertex = _v[0] == vertex ? 0 : 1; }
virtual int getNumEdges(){ return 1; }
virtual MEdge getEdge(int num){ return MEdge(_v[0], _v[1]); }
virtual int getNumEdgesRep(){ return 1; }
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 55ed9c2c69..3098159e06 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -50,10 +50,22 @@ class MPoint : public MElement {
{
s[0][0] = s[0][1] = s[0][2] = 0.;
}
+
+ virtual const polynomialBasis* getFunctionSpace(int o) const { return &polynomialBases::find(MSH_PNT); }
virtual bool isInside(double u, double v, double w)
{
return true;
}
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
+ {
+ static IntPt GQL[1];
+ GQL[0].pt[0] = 0;
+ GQL[0].pt[1] = 0;
+ GQL[0].pt[2] = 0;
+ GQL[0].weight = 1;
+ *npts = 1;
+ *pts = GQL;
+ }
};
#endif
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 46f7c9e38e..69c28525df 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -21,10 +21,13 @@ static fullMatrix<double> generate1DMonomials(int order)
static fullMatrix<double> generate1DPoints(int order)
{
fullMatrix<double> line(order + 1, 1);
- line(0, 0) = -1.;
- line(1, 0) = 1.;
- double dd = 2. / order;
- for (int i = 2; i < order + 1; i++) line(i, 0) = -1. + dd * (i - 1);
+ line(0,0) = 0;
+ if (order > 0) {
+ line(0, 0) = -1.;
+ line(1, 0) = 1.;
+ double dd = 2. / order;
+ for (int i = 2; i < order + 1; i++) line(i, 0) = -1. + dd * (i - 1);
+ }
return line;
}
@@ -762,6 +765,14 @@ static void generate2dEdgeClosure(polynomialBasis::clCont &closure, int order)
}
}
+
+static void generate1dVertexClosure(polynomialBasis::clCont &closure)
+{
+ closure.clear();
+ closure.resize(2);
+ closure[0].push_back(0);
+ closure[1].push_back(1);
+ }
std::map<int, polynomialBasis> polynomialBases::fs;
const polynomialBasis &polynomialBases::find(int tag)
@@ -779,22 +790,27 @@ const polynomialBasis &polynomialBases::find(int tag)
case MSH_LIN_2 :
F.monomials = generate1DMonomials(1);
F.points = generate1DPoints(1);
+ generate1dVertexClosure(F.vertexClosure);
break;
case MSH_LIN_3 :
F.monomials = generate1DMonomials(2);
F.points = generate1DPoints(2);
+ generate1dVertexClosure(F.vertexClosure);
break;
case MSH_LIN_4:
F.monomials = generate1DMonomials(3);
F.points = generate1DPoints(3);
+ generate1dVertexClosure(F.vertexClosure);
break;
case MSH_LIN_5:
F.monomials = generate1DMonomials(4);
F.points = generate1DPoints(4);
+ generate1dVertexClosure(F.vertexClosure);
break;
case MSH_LIN_6:
F.monomials = generate1DMonomials(5);
F.points = generate1DPoints(5);
+ generate1dVertexClosure(F.vertexClosure);
break;
case MSH_TRI_3 :
F.monomials = generatePascalTriangle(1);
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index 01b8f846c3..9abcf5549e 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -18,6 +18,7 @@ struct polynomialBasis
typedef std::vector<std::vector<int> > clCont;
clCont faceClosure;
clCont edgeClosure;
+ clCont vertexClosure;
fullMatrix<double> points;
fullMatrix<double> monomials;
fullMatrix<double> coefficients;
@@ -31,6 +32,10 @@ struct polynomialBasis
{
return edgeClosure[iSign == 1 ? iEdge : 3 + iEdge];
}
+ inline const std::vector<int> &getVertexClosure(int iVertex) const
+ {
+ return vertexClosure[iVertex];
+ }
inline void evaluateMonomials(double u, double v, double w, double p[]) const
{
for (int j = 0; j < monomials.size1(); j++) {
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index 9a1345a569..9d24b82703 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -364,10 +364,9 @@ void dgAlgorithm::buildGroups(GModel *model, int dim, int order,
switch(dim) {
case 1 : {
std::map<const std::string, std::set<MVertex*> >::iterator mapIt;
- /*for(mapIt=boundaryVertices.begin(); mapIt!=boundaryVertices.end(); mapIt++) {
+ for(mapIt=boundaryVertices.begin(); mapIt!=boundaryVertices.end(); mapIt++) {
bGroups.push_back(new dgGroupOfFaces(*eGroups[0],mapIt->first,order,mapIt->second));
- }*/
- throw;
+ }
break;
}
case 2 : {
diff --git a/Solver/dgGroupOfElements.cpp b/Solver/dgGroupOfElements.cpp
index 7694e1c065..315b037103 100644
--- a/Solver/dgGroupOfElements.cpp
+++ b/Solver/dgGroupOfElements.cpp
@@ -4,6 +4,7 @@
#include "Numeric.h"
#include "MTriangle.h"
#include "MLine.h"
+#include "MPoint.h"
#include "GModel.h"
static fullMatrix<double> * dgGetIntegrationRule (MElement *e, int p){
@@ -68,16 +69,15 @@ dgGroupOfElements::dgGroupOfElements(const std::vector<MElement*> &e, int polyOr
e->getJacobian ((*_integration)(j,0), (*_integration)(j,1), (*_integration)(j,2), jac);
const double weight = (*_integration)(j,3);
detjac=inv3x3(jac,ijac);
- (*_mapping)(i,10*j + 0) = ijac[0][0];
- (*_mapping)(i,10*j + 1) = ijac[0][1];
- (*_mapping)(i,10*j + 2) = ijac[0][2];
- (*_mapping)(i,10*j + 3) = ijac[1][0];
- (*_mapping)(i,10*j + 4) = ijac[1][1];
- (*_mapping)(i,10*j + 5) = ijac[1][2];
- (*_mapping)(i,10*j + 6) = ijac[2][0];
- (*_mapping)(i,10*j + 7) = ijac[2][1];
- (*_mapping)(i,10*j + 8) = ijac[2][2];
- (*_mapping)(i,10*j + 9) = detjac;
+ (*_mapping)(i,10*j + 0) = ijac[0][0];
+ (*_mapping)(i,10*j + 1) = ijac[0][1];
+ (*_mapping)(i,10*j + 2) = ijac[0][2];
+ (*_mapping)(i,10*j + 3) = ijac[1][0];
+ (*_mapping)(i,10*j + 4) = ijac[1][1];
+ (*_mapping)(i,10*j + 5) = ijac[1][2];
+ (*_mapping)(i,10*j + 6) = ijac[2][0];
+ (*_mapping)(i,10*j + 7) = ijac[2][1];
+ (*_mapping)(i,10*j + 8) = ijac[2][2];
for (int k=0;k<_fs.coefficients.size1();k++){
for (int l=0;l<_fs.coefficients.size1();l++) {
imass(k,l) += f[k]*f[l]*weight*detjac;
@@ -129,8 +129,8 @@ void dgGroupOfFaces::addFace(const MFace &topoFace, int iElLeft, int iElRight){
_closuresLeft.push_back(&(_fsLeft->getFaceClosure(ithFace, sign, rot)));
const std::vector<int> & geomClosure = elLeft.getFunctionSpace()->getFaceClosure(ithFace, sign, rot);
if(iElRight>=0){
- MElement &elRight = *_groupRight.getElement(iElRight);
_right.push_back(iElRight);
+ MElement &elRight = *_groupRight.getElement(iElRight);
elRight.getFaceInfo (topoFace, ithFace, sign, rot);
_closuresRight.push_back(&(_fsRight->getFaceClosure(ithFace, sign, rot)));
}
@@ -185,6 +185,21 @@ void dgGroupOfFaces::addEdge(const MEdge &topoEdge, int iElLeft, int iElRight){
}
}
+void dgGroupOfFaces::addVertex(MVertex *topoVertex, int iElLeft, int iElRight){
+ _left.push_back(iElLeft);
+ MElement &elLeft = *_groupLeft.getElement(iElLeft);
+ int ithVertex;
+ elLeft.getVertexInfo (topoVertex, ithVertex);
+ _closuresLeft.push_back(&_fsLeft->getVertexClosure(ithVertex));
+ if(iElRight>=0){
+ _right.push_back(iElRight);
+ MElement &elRight = *_groupRight.getElement(iElRight);
+ elRight.getVertexInfo (topoVertex, ithVertex);
+ _closuresRight.push_back(&_fsRight->getVertexClosure(ithVertex));
+ }
+ _faces.push_back(new MPoint(topoVertex) );
+}
+
void dgGroupOfFaces::init(int pOrder) {
_fsFace = _faces[0]->getFunctionSpace (pOrder);
_integration = dgGetIntegrationRule (_faces[0],pOrder);
@@ -283,6 +298,26 @@ dgGroupOfFaces::~dgGroupOfFaces()
delete _dPsiLeftDxOnQP;
delete _dPsiRightDxOnQP;
}
+
+dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices):
+ _groupLeft(elGroup),_groupRight(elGroup)
+{
+ _boundaryTag=boundaryTag;
+ if(boundaryTag=="")
+ throw;
+ _fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
+ _fsRight=NULL;
+ for(int i=0; i<elGroup.getNbElements(); i++){
+ MElement &el = *elGroup.getElement(i);
+ for (int j=0; j<el.getNumVertices(); j++){
+ MVertex* vertex = el.getVertex(j);
+ if(boundaryVertices.find(vertex)!=boundaryVertices.end())
+ addVertex(vertex,i,-1);
+ }
+ }
+ init(pOrder);
+}
+
dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MEdge,Less_Edge> &boundaryEdges):
_groupLeft(elGroup),_groupRight(elGroup)
{
@@ -301,6 +336,7 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string bo
}
init(pOrder);
}
+
dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MFace,Less_Face> &boundaryFaces):
_groupLeft(elGroup),_groupRight(elGroup)
{
@@ -326,6 +362,21 @@ dgGroupOfFaces::dgGroupOfFaces (const dgGroupOfElements &elGroup, int pOrder):
_fsLeft=_groupLeft.getElement(0)->getFunctionSpace(pOrder);
_fsRight=_groupRight.getElement(0)->getFunctionSpace(pOrder);
switch (_groupLeft.getElement(0)->getDim()) {
+ case 1 : {
+ std::map<MVertex*,int> vertexMap;
+ for(int i=0; i<elGroup.getNbElements(); i++){
+ MElement &el = *elGroup.getElement(i);
+ for (int j=0; j<el.getNumVertices(); j++){
+ MVertex* vertex = el.getVertex(j);
+ if(vertexMap.find(vertex) == vertexMap.end()){
+ vertexMap[vertex] = i;
+ }else{
+ addVertex(vertex,vertexMap[vertex],i);
+ }
+ }
+ }
+ break;
+ }
case 2 : {
std::map<MEdge,int,Less_Edge> edgeMap;
for(int i=0; i<elGroup.getNbElements(); i++){
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 3482bbb4da..2bd2e554aa 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -90,6 +90,7 @@ class dgGroupOfFaces {
const dgGroupOfElements &_groupLeft,&_groupRight;
void addFace(const MFace &topoFace, int iElLeft, int iElRight);
void addEdge(const MEdge &topoEdge, int iElLeft, int iElRight);
+ void addVertex(MVertex *topoVertex, int iElLeft, int iElRight);
// Two polynomialBases for left and right elements
// the group has always the same types for left and right
const polynomialBasis *_fsLeft,*_fsRight, *_fsFace;
@@ -132,6 +133,7 @@ public:
const std::vector<int> * getClosureRight(int iFace) const{ return _closuresRight[iFace];}
inline fullMatrix<double> &getNormals () const {return *_normals;}
dgGroupOfFaces (const dgGroupOfElements &elements,int pOrder);
+ dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MVertex*> &boundaryVertices);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MEdge,Less_Edge> &boundaryEdges);
dgGroupOfFaces (const dgGroupOfElements &elGroup, std::string boundaryTag, int pOrder,std::set<MFace,Less_Face> &boundaryFaces);
virtual ~dgGroupOfFaces ();
--
GitLab