From ec0b07dbf036d67adbab21656c27ebde28e35c7f Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 20 Jan 2010 19:21:59 +0000
Subject: [PATCH] smoothing high order works with quads !!
---
Common/LuaBindings.cpp | 1 +
Geo/GModel.cpp | 4 ++
Geo/MElement.cpp | 16 +++++
Geo/MElement.h | 2 +
Geo/MQuadrangle.cpp | 26 ++++++++
Geo/MQuadrangle.h | 2 +
Geo/MTriangle.cpp | 1 +
Geo/MVertex.cpp | 2 +-
Mesh/HighOrder.cpp | 15 ++++-
Mesh/highOrderSmoother.cpp | 2 +-
Mesh/qualityMeasures.cpp | 92 ++++++++++++++++++++++++++++-
Mesh/qualityMeasures.h | 3 +
Solver/TESTCASES/CylinderEddies.lua | 4 +-
Solver/TESTCASES/cyl2.geo | 8 +--
Solver/dgAlgorithm.cpp | 12 ++--
Solver/dgGroupOfElements.h | 1 +
Solver/dgSystemOfEquations.h | 3 +
Solver/function.cpp | 38 ++++++++++++
Solver/function.h | 10 ++++
19 files changed, 225 insertions(+), 17 deletions(-)
diff --git a/Common/LuaBindings.cpp b/Common/LuaBindings.cpp
index af0d3cf8fe..1b98e8ca00 100644
--- a/Common/LuaBindings.cpp
+++ b/Common/LuaBindings.cpp
@@ -179,6 +179,7 @@ binding::binding(){
functionLua::registerBindings(this);
function::registerDefaultFunctions();
MVertex::registerBindings(this);
+ MElement::registerBindings(this);
}
binding *binding::_instance=NULL;
#endif
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 74d614020c..0bd0b84f3d 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1401,5 +1401,9 @@ void GModel::registerBindings(binding *b)
cm = cb->addMethod("mesh", &GModel::mesh);
cm = cb->addMethod("load", &GModel::load);
cm = cb->addMethod("save", &GModel::save);
+ cm = cb->addMethod("getNumMeshElements",(int (GModel::*)()) &GModel::getNumMeshElements);
+ cm = cb->addMethod("getMeshElementByTag",&GModel::getMeshElementByTag);
+ cm = cb->addMethod("getNumMeshVertices",&GModel::getNumMeshVertices);
+ cm = cb->addMethod("getMeshVertexByTag",&GModel::getMeshVertexByTag);
cm = cb->setConstructor<GModel>();
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index b29d77a43b..be20d7389b 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -801,3 +801,19 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
default: return 0;
}
}
+
+#include "Bindings.h"
+
+void MElement::registerBindings(binding *b)
+{
+ classBinding *cb = b->addClass<MElement>("MElement");
+ methodBinding *cm;
+ cm = cb->addMethod("getNum",&MElement::getNum);
+ // here we specify the cast because there are 2 MVertex::x function
+ cm = cb->addMethod("getNumVertices", &MElement::getNumVertices);
+ cm = cb->addMethod("getVertex", &MElement::getVertex);
+ cm = cb->addMethod("getType", &MElement::getType);
+ cm = cb->addMethod("getPartition", &MElement::getPartition);
+ cm = cb->addMethod("getPolynomialOrder", &MElement::getPolynomialOrder);
+ cm = cb->addMethod("getDim", &MElement::getDim);
+}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 71df81d183..fbd38e5616 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -18,6 +18,7 @@
#include "Gauss.h"
class GFace;
+class binding;
// A mesh element.
class MElement
@@ -266,6 +267,7 @@ class MElement
// return the number of vertices, as well as the element name if
// 'name' != 0
static int getInfoMSH(const int typeMSH, const char **const name=0);
+ static void registerBindings(binding *b);
};
class MElementLessThanLexicographic{
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index aec951b0c4..955a2507d0 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -7,6 +7,11 @@
#include "GaussLegendre1D.h"
#include "Context.h"
#include "qualityMeasures.h"
+
+#if defined(HAVE_MESH)
+#include "qualityMeasures.h"
+#endif
+
#define SQU(a) ((a)*(a))
const polynomialBasis* MQuadrangle::getFunctionSpace(int o) const
{
@@ -38,6 +43,8 @@ int MQuadrangleN::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges
int MQuadrangle8::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
int MQuadrangle9::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
+
+
static void _myGetEdgeRep(MQuadrangle *q, int num, double *x, double *y, double *z,
SVector3 *n, int numSubEdges)
{
@@ -159,3 +166,22 @@ void MQuadrangle::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
*pts = getGQQPts(pOrder);
}
+double MQuadrangle::distoShapeMeasure()
+{
+#if defined(HAVE_MESH)
+ // return qmTriangleAngles(this);
+ return qmDistorsionOfMapping(this);
+#else
+ return 0.;
+#endif
+}
+
+
+double MQuadrangle::angleShapeMeasure()
+{
+#if defined(HAVE_MESH)
+ return qmQuadrangleAngles(this);
+#else
+ return 1.;
+#endif
+}
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index a4c274f6d0..ecfe3d254e 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -129,6 +129,8 @@ class MQuadrangle : public MElement {
return true;
}
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
+ virtual double angleShapeMeasure();
+ virtual double distoShapeMeasure();
private:
int edges_quad(const int edge, const int vert) const
{
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index ab7171141e..8b4a4a9a40 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -27,6 +27,7 @@ SPoint3 MTriangle::circumcenter()
double MTriangle::distoShapeMeasure()
{
#if defined(HAVE_MESH)
+ //return qmTriangleAngles(this);
return qmDistorsionOfMapping(this);
#else
return 0.;
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 0f56dfdd79..c0e19b44e5 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -407,7 +407,7 @@ void MVertex::registerBindings(binding *b)
cm = cb->addMethod("getNum",&MVertex::getNum);
// here we specify the cast because there are 2 MVertex::x function
cm = cb->addMethod("x", (double (MVertex::*)() const) &MVertex::x);
- cm = cb->addMethod("x2", (double& (MVertex::*)()) &MVertex::x);
+ // cm = cb->addMethod("x2", (double& (MVertex::*)()) &MVertex::x);
cm = cb->addMethod("y", (double (MVertex::*)() const) &MVertex::y);
cm = cb->addMethod("z", (double (MVertex::*)() const) &MVertex::z);
cm = cb->setConstructor<MVertex,double,double,double>();
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index b36767ecd9..5ab8d271dd 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -942,12 +942,23 @@ static void checkHighOrderTriangles(const char* cc, GModel *m,
if (disto < 0) bad.push_back(t);
else if (disto < 0.2) nbfair++;
}
+ for(unsigned int i = 0; i < (*it)->quadrangles.size(); i++){
+ MQuadrangle *t = (*it)->quadrangles[i];
+ double disto_ = t->distoShapeMeasure();
+ double gamma_ = t->gammaShapeMeasure();
+ double disto = disto_;
+ minJGlob = std::min(minJGlob, disto);
+ minGGlob = std::min(minGGlob, gamma_);
+ avg += disto; count++;
+ if (disto < 0) bad.push_back(t);
+ else if (disto < 0.2) nbfair++;
+ }
}
if (minJGlob > 0)
- Msg::Info("%s : Worst Triangle Smoothness %g Gamma %g NbFair = %d",
+ Msg::Info("%s : Worst Face Smoothness %g Gamma %g NbFair = %d",
cc, minJGlob, minGGlob,nbfair );
else
- Msg::Warning("%s : Worst Triangle Smoothness %g (%d negative jacobians) "
+ Msg::Warning("%s : Worst Face Smoothness %g (%d negative jacobians) "
"Worst Gamma %g Avg Smoothness %g", cc, minJGlob, bad.size(),
minGGlob, avg / (count ? count : 1));
}
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index a885beeae7..683f3e5831 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -369,7 +369,7 @@ void highOrderSmoother::optimize(GFace * gf,
smooth(gf, true);
//int nbSwap =
- swapHighOrderTriangles(gf,edgeVertices,faceVertices,this);
+ // swapHighOrderTriangles(gf,edgeVertices,faceVertices,this);
// smooth(gf,true);
// smooth(gf,true);
// smooth(gf,true);
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 330cb4f852..bfa40d89bc 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -7,6 +7,7 @@
#include "BDS.h"
#include "MVertex.h"
#include "MTriangle.h"
+#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "Numeric.h"
#include "polynomialBasis.h"
@@ -177,7 +178,7 @@ double qmTet(const double &x1, const double &y1, const double &z1,
}
}
-double mesh_functional_distorsion(MTriangle *t, double u, double v)
+double mesh_functional_distorsion(MElement *t, double u, double v)
{
// compute uncurved element jacobian d_u x and d_v x
double mat[3][3];
@@ -247,6 +248,35 @@ double qmDistorsionOfMapping (MTriangle *e)
return dmin;
}
+double qmDistorsionOfMapping (MQuadrangle *e)
+{
+ // return 1.0;
+ if (e->getPolynomialOrder() == 1) return 1.0;
+ // if (e->getPolynomialOrder() == 2) return mesh_functional_distorsion_p2(e);
+
+ IntPt *pts;
+ int npts;
+ e->getIntegrationPoints(e->getPolynomialOrder(),&npts, &pts);
+ double dmin;
+ for (int i = 0 ; i < npts; i++){
+ const double u = pts[i].pt[0];
+ const double v = pts[i].pt[1];
+ const double di = mesh_functional_distorsion (e, u, v);
+ dmin = (i == 0)? di : std::min(dmin, di);
+ }
+ const fullMatrix<double>& points = e->getFunctionSpace()->points;
+
+ for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
+ const double u = points(i, 0);
+ const double v = points(i, 1);
+ const double di = mesh_functional_distorsion (e, u, v);
+ dmin = std::min(dmin, di);
+ }
+ // printf("%g\n",dmin);
+ return dmin;
+}
+
+
static double mesh_functional_distorsion(MTetrahedron *t, double u, double v, double w)
{
// compute uncurved element jacobian d_u x and d_v x
@@ -351,3 +381,63 @@ double qmTriangleAngles (MTriangle *e) {
return worst_quality;
}
+
+
+double qmQuadrangleAngles (MQuadrangle *e) {
+ double a = 100;
+ double worst_quality = std::numeric_limits<double>::max();
+ double mat[3][3];
+ double mat2[3][3];
+ double den = atan(a*(M_PI/4)) + atan(a*(2*M_PI/4 - (M_PI/4)));
+
+ // This matrix is used to "rotate" the triangle to get each vertex
+ // as the "origin" of the mapping in turn
+ double rot[3][3];
+ rot[0][0]=-1; rot[0][1]=1; rot[0][2]=0;
+ rot[1][0]=-1; rot[1][1]=0; rot[1][2]=0;
+ rot[2][0]= 0; rot[2][1]=0; rot[2][2]=1;
+ double tmp[3][3];
+
+ const double u[9] = {-1,-1, 1, 1, 0,0,1,-1,0};
+ const double v[9] = {-1, 1, 1,-1, -1,1,0,0,0};
+
+ for (int i = 0; i < 9; i++) {
+
+ e->getJacobian(u[i], v[i], 0, mat);
+ e->getPrimaryJacobian(u[i],v[i],0,mat2);
+ // for (int j = 0; j < i; j++) {
+ //matmat(rot,mat,tmp);
+ // memcpy(mat, tmp, sizeof(mat));
+ // }
+ //get angle
+ double v1[3] = {mat[0][0], mat[0][1], mat[0][2] };
+ double v2[3] = {mat[1][0], mat[1][1], mat[1][2] };
+ double v3[3] = {mat2[0][0], mat2[0][1], mat2[0][2] };
+ double v4[3] = {mat2[1][0], mat2[1][1], mat2[1][2] };
+ norme(v1);
+ norme(v2);
+ norme(v3);
+ norme(v4);
+ double v12[3], v34[3];
+ prodve(v1,v2,v12);
+ prodve(v3,v4,v34);
+ norme(v12);
+ norme(v34);
+ double orientation;
+ prosca(v12,v34,&orientation);
+
+ // If the if the triangle is "flipped" it's no good
+ // if (orientation < 0)
+ // return -std::numeric_limits<double>::max();
+
+ double c;
+ prosca(v1,v2,&c);
+ printf("%g %g\n",c,acos(c)*180/M_PI);
+ double x = fabs(acos(c))-M_PI/2;
+ double quality = (atan(a*(x+M_PI/4)) + atan(a*(2*M_PI/4 - (x+M_PI/4))))/den;
+ worst_quality = std::min(worst_quality, quality);
+ }
+
+ return worst_quality;
+}
+
diff --git a/Mesh/qualityMeasures.h b/Mesh/qualityMeasures.h
index cacf55add0..16be388eb1 100644
--- a/Mesh/qualityMeasures.h
+++ b/Mesh/qualityMeasures.h
@@ -10,16 +10,19 @@ class BDS_Point;
class BDS_Face;
class MVertex;
class MTriangle;
+class MQuadrangle;
class MTetrahedron;
class MElement;
enum qualityMeasure4Triangle {QMTRI_RHO, QMTRI_COND};
enum qualityMeasure4Tet{QMTET_1, QMTET_2, QMTET_3, QMTET_ONE, QMTET_COND};
+double qmDistorsionOfMapping (MQuadrangle *e);
double qmDistorsionOfMapping (MTriangle *e);
double qmDistorsionOfMapping (MTetrahedron *e);
double qmTriangleAngles (MTriangle *e);
+double qmQuadrangleAngles (MQuadrangle *e);
double qmTriangle(MTriangle *f, const qualityMeasure4Triangle &cr);
double qmTriangle(BDS_Face *f, const qualityMeasure4Triangle &cr);
diff --git a/Solver/TESTCASES/CylinderEddies.lua b/Solver/TESTCASES/CylinderEddies.lua
index 42b5d087ce..0a7d9cf3c5 100644
--- a/Solver/TESTCASES/CylinderEddies.lua
+++ b/Solver/TESTCASES/CylinderEddies.lua
@@ -20,7 +20,7 @@ end
Example of a lua program driving the DG code
--]]
-order = 3
+order = 5
FS = functionLua(4, 'free_stream', {'XYZ'}):getName()
@@ -62,7 +62,7 @@ CFL = 20.1;
dt = CFL * DG:computeInvSpectralRadius();
print('DT = ',dt)
T = 0;
-for i=1,100000 do
+for i=1,1000000 do
dt = CFL * DG:computeInvSpectralRadius();
norm = DG:RK44(dt)
T = T + dt
diff --git a/Solver/TESTCASES/cyl2.geo b/Solver/TESTCASES/cyl2.geo
index 08735cab6f..1e8ee12b1e 100644
--- a/Solver/TESTCASES/cyl2.geo
+++ b/Solver/TESTCASES/cyl2.geo
@@ -1,11 +1,11 @@
-a = .1;
-b = 1;
+a = .4;
+b = 2;
radius = 0.5;
radiusBoundaryLayer = 0.7;
Point(1) = {0, 0, 0, a};
-nlayers=3;
-np=5;
+nlayers=2;
+np=3;
np2=2;
diff --git a/Solver/dgAlgorithm.cpp b/Solver/dgAlgorithm.cpp
index fd3ad00617..8035d30533 100644
--- a/Solver/dgAlgorithm.cpp
+++ b/Solver/dgAlgorithm.cpp
@@ -117,12 +117,12 @@ void dgAlgorithm::residualVolume ( //dofManager &dof, // the DOF manager (maybe
}
void dgAlgorithm::residualInterface ( //dofManager &dof, // the DOF manager (maybe useless here)
- const dgConservationLaw &claw, // the conservation law
- dgGroupOfFaces &group,
- const fullMatrix<double> &solution, // solution !! at faces nodes
- fullMatrix<double> &solutionLeft,
- fullMatrix<double> &solutionRight,
- fullMatrix<double> &residual // residual !! at faces nodes
+ const dgConservationLaw &claw, // the conservation law
+ dgGroupOfFaces &group,
+ const fullMatrix<double> &solution, // solution !! at faces nodes
+ fullMatrix<double> &solutionLeft,
+ fullMatrix<double> &solutionRight,
+ fullMatrix<double> &residual // residual !! at faces nodes
)
{
// A) global operations before entering the loop over the faces
diff --git a/Solver/dgGroupOfElements.h b/Solver/dgGroupOfElements.h
index 55f742daa4..0b9be6a2d1 100644
--- a/Solver/dgGroupOfElements.h
+++ b/Solver/dgGroupOfElements.h
@@ -82,6 +82,7 @@ public:
inline int getDimUVW () const {return _dimUVW;}
inline int getDimXYZ () const {return _dimXYZ;}
inline MElement* getElement (int iElement) const {return _elements[iElement];}
+ inline const polynomialBasis & getFunctionSpace () const {return _fs;}
inline const fullMatrix<double> & getIntegrationPointsMatrix () const {return *_integration;}
inline const fullMatrix<double> & getCollocationMatrix () const {return *_collocation;}
inline const fullMatrix<double> & getFluxRedistributionMatrix (int i) const {return *_redistributionFluxes[i];}
diff --git a/Solver/dgSystemOfEquations.h b/Solver/dgSystemOfEquations.h
index e8d89dff59..fdb887dc93 100644
--- a/Solver/dgSystemOfEquations.h
+++ b/Solver/dgSystemOfEquations.h
@@ -36,6 +36,9 @@ class dgSystemOfEquations {
dgSystemOfEquations(const dgSystemOfEquations &) {}
double computeTimeStepMethodOfLines () const;
public:
+ const dgConservationLaw * getLaw() const {return _claw;}
+ const GModel * getModel() const {return _gm;}
+ std::pair<dgGroupOfElements*,int> getElementPosition (MElement *);
void setOrder (int order); // set the polynomial order
void setConservationLaw (dgConservationLaw *law); // set the conservationLaw
dgSystemOfEquations(GModel *_gm);
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 00a19a8b9e..998d8358d5 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -143,6 +143,44 @@ functionConstant::functionConstant(const fullMatrix<double> *source){
function::add(_name,this);
}
+// function that enables to interpolate a DG solution using
+// geometrical search in a mesh
+/*
+class functionSystemOfEquations::data : public dataCacheDouble {
+ dgSystemOfEquations *_sys;
+ dataCacheDouble &xyz;
+public:
+ data(dataCacheMap &m, dgSystemOfEquations *sys) :
+ _sys(sys),
+ xyz(cacheMap.get("Solution",this)),
+ dataCacheDouble(m.getNbEvaluationPoints(), sys->getLaw()->getNbFields())
+ {
+ }
+ void _eval() {
+ int nP =xyz().size1();
+ if(_value.size1() != nP)
+ _value = fullMatrix<double>(nP,sys->getLaw()->getNbFields());
+ _value.setAll(0.0);
+ double fs[256];
+ for (int i=0;i<_value.size1();i++){
+ const double x = xyz(i,0);
+ const double y = xyz(i,1);
+ const double z = xyz(i,2);
+ MElement *e = _sys->getModel()->getMeshElementByCoord(SPoint3(x,y,z));
+ std::pair<dgGroupOfElements*,int> location = _sys->getElementPosition(e);
+ double U[3],X[3]={xyz(i,0),xyz(i,1),xyz(i,2)};
+ e->xyz2uvw (X,U);
+ location.first->getFunctionSpace().f(U[0],U[1],U[2],fs);
+ }
+ }
+ dataCacheDouble *newDataCache(dataCacheMap *m)
+ {
+ return new data(this,_sys);
+ }
+};
+*/
+
+
#include "Bindings.h"
void function::registerDefaultFunctions()
diff --git a/Solver/function.h b/Solver/function.h
index 51916a46d1..9a81cbfa7d 100644
--- a/Solver/function.h
+++ b/Solver/function.h
@@ -8,6 +8,7 @@
class dataCacheMap;
class MElement;
class binding;
+class dgSystemOfEquations;
// those classes manage complex function dependencies and keep their values in cache so that they are not recomputed when it is not necessary. To do this, we use three classes : function, dataCache and dataCacheMap. The workflow is :
//
@@ -193,4 +194,13 @@ class functionConstant : public function {
dataCacheDouble *newDataCache(dataCacheMap *m);
functionConstant(const fullMatrix<double> *source);
};
+
+class functionSystemOfEquations : public function {
+ dgSystemOfEquations *_sys;
+public:
+ class data ;
+ dataCacheDouble *newDataCache(dataCacheMap *m);
+ functionSystemOfEquations(dgSystemOfEquations *sys) : _sys(sys) {}
+};
+
#endif
--
GitLab