From bf66c842fb8a492b13e7709f48d68b74235bb6c2 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 5 Apr 2011 10:39:18 +0000
Subject: [PATCH] mesh matching enhanced high order smoother class added : high
order meshes are guaranteed to be valid : some elements are "un-curved" in
order to ensure good value of the distorsion mapping. P2 triangle quality
estimate is robust
---
Common/gmshpy.i | 2 +
Fltk/classificationEditor.cpp | 12 +-
Fltk/classificationEditor.h | 2 +-
Geo/GEdge.cpp | 97 +++--
Geo/GModelFactory.cpp | 3 +-
Geo/GeomMeshMatcher.cpp | 169 +++++++++
Geo/GeomMeshMatcher.h | 2 +
Geo/MElement.h | 10 +
Mesh/CMakeLists.txt | 1 +
Mesh/HighOrder.cpp | 24 +-
Mesh/highOrderSmoother.cpp | 77 +++-
Mesh/highOrderSmoother.h | 1 -
Mesh/highOrderTools.cpp | 668 ++++++++++++++++++++++++++++++++++
Mesh/highOrderTools.h | 92 +++++
Mesh/meshGFaceOptimize.cpp | 6 +
Mesh/meshGFaceOptimize.h | 2 +
Mesh/qualityMeasures.cpp | 18 +-
Numeric/DivideAndConquer.cpp | 55 ++-
Numeric/DivideAndConquer.h | 18 +-
Solver/elasticityTerm.cpp | 71 +++-
20 files changed, 1253 insertions(+), 77 deletions(-)
create mode 100644 Mesh/highOrderTools.cpp
create mode 100644 Mesh/highOrderTools.h
diff --git a/Common/gmshpy.i b/Common/gmshpy.i
index 452afd9951..0c87dfb7f1 100644
--- a/Common/gmshpy.i
+++ b/Common/gmshpy.i
@@ -6,6 +6,7 @@
%{
#include "GmshConfig.h"
#include "GModel.h"
+ #include "highOrderTools.h"
#include "DefaultOptions.h"
#include "fullMatrix.h"
#include "function.h"
@@ -92,6 +93,7 @@ namespace std {
%include "dofManager.h"
%template(dofManagerDouble) dofManager<double>;
%include "GModel.h"
+%include "highOrderTools.h"
%include "function.h"
%include "linearSystem.h"
%template(linearSystemDouble) linearSystem<double>;
diff --git a/Fltk/classificationEditor.cpp b/Fltk/classificationEditor.cpp
index 7c76c9c505..b9203faf09 100644
--- a/Fltk/classificationEditor.cpp
+++ b/Fltk/classificationEditor.cpp
@@ -15,6 +15,7 @@
#include "Context.h"
#include "GmshMessage.h"
#include "MLine.h"
+#include "MQuadrangle.h"
#include "meshGFaceDelaunayInsertion.h"
#include "discreteEdge.h"
#include "discreteFace.h"
@@ -104,9 +105,12 @@ static void select_elements_cb(Fl_Widget *w, void *data)
if(all){
for(GModel::fiter it = GModel::current()->firstFace();
- it != GModel::current()->lastFace(); ++it)
+ it != GModel::current()->lastFace(); ++it){
e->elements.insert(e->elements.end(), (*it)->triangles.begin(),
(*it)->triangles.end());
+ e->elements.insert(e->elements.end(), (*it)->quadrangles.begin(),
+ (*it)->quadrangles.end());
+ }
}
else{
CTX::instance()->pickElements = 1;
@@ -120,9 +124,9 @@ static void select_elements_cb(Fl_Widget *w, void *data)
if(ib == 'l') {
for(unsigned int i = 0; i < FlGui::instance()->selectedElements.size(); i++){
MElement *me = FlGui::instance()->selectedElements[i];
- if(me->getType() == TYPE_TRI && me->getVisibility() != 2){
+ if(me->getDim() == 2 && me->getVisibility() != 2){
me->setVisibility(2);
- e->elements.push_back((MTriangle*)me);
+ e->elements.push_back(me);
}
}
}
@@ -148,7 +152,7 @@ static void select_elements_cb(Fl_Widget *w, void *data)
}
e2t_cont adj;
- buildEdgeToTriangle(e->elements, adj);
+ buildEdgeToElements(e->elements, adj);
buildListOfEdgeAngle(adj, e->edges_detected, e->edges_lonly);
ElementsSelectedMode(e);
update_edges_cb(0, data);
diff --git a/Fltk/classificationEditor.h b/Fltk/classificationEditor.h
index c265da22c8..d7eee0113f 100644
--- a/Fltk/classificationEditor.h
+++ b/Fltk/classificationEditor.h
@@ -31,7 +31,7 @@
class classificationEditor {
public:
- std::vector<MTriangle*> elements;
+ std::vector<MElement*> elements;
std::set<GFace*> faces;
Fl_Window *window;
Fl_Button *buttons[10];
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index adf7d5431c..7e0a86b73b 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -259,38 +259,83 @@ double GEdge::length(const double &u0, const double &u1, const int nbQuadPoints)
return L;
}
+/*
+ consider a curve x(t) and a point y
+
+ use a golden section algorithm that minimizes
+
+ min_t \|x(t)-y\|
+
+*/
+
+const double GOLDEN = (1. + sqrt(5)) / 2.;
+const double GOLDEN2 = 2 - GOLDEN;
+
+// x1 and x3 are the current bounds; the minimum is between them.
+// x2 is the center point, which is closer to x1 than to x3
+
+double goldenSectionSearch(const GEdge *ge, const SPoint3 &q, double x1, double x2, double x3, double tau)
+{
+ // Create a new possible center in the area between x2 and x3, closer to x2
+ double x4 = x2 + GOLDEN2 * (x3 - x2);
+
+ // Evaluate termination criterion
+ if (fabs(x3 - x1) < tau * (fabs(x2) + fabs(x4)))
+ return (x3 + x1) / 2;
+
+ const SVector3 dp4 = q - ge->position(x4);
+ const SVector3 dp2 = q - ge->position(x2);
+
+ const double d4 = dp4.norm();
+ const double d2 = dp2.norm();
+
+ if (d4 < d2)
+ return goldenSectionSearch(ge, q, x2, x4, x3, tau);
+ else
+ return goldenSectionSearch(ge,q , x4, x2, x1, tau);
+}
+
GPoint GEdge::closestPoint(const SPoint3 &q, double &t) const
{
+
+ // printf("looking for closest point in curve %d to point %g %g\n",tag(),q.x(),q.y());
+
+ const int nbSamples = 100;
+
double tolerance = 1.e-12;
Range<double> interval = parBounds(0);
double tMin = std::min(interval.high(), interval.low());
double tMax = std::max(interval.high(), interval.low());
- double relax = 1.;
- double dt, dt0, t0;
- int nb = 10;
-
- t = (tMin + tMax) * 0.5;
-
- while (relax > 0.1) {
- int i = 0;
- t = 0.5 * (tMin + tMax);
- dt0 = tMax - tMin;
- dt = dt0;
- while (dt > tolerance * dt0 && i++ < nb) {
- t0 = t;
- dt0 = dt;
- SVector3 dp = q - position(t);
- SVector3 derP = firstDer(t);
- double b = dot(derP, derP);
- double c = dot(derP, dp);
- dt = c / b;
- t = std::max(tMin, std::min(tMax, t0 + relax * dt));
- dt = fabs(t - t0);
+
+ double DMIN = 1.e22;
+ double topt;
+ const double DT = (tMax-tMin)/(nbSamples-1) ;
+ for (int i=0;i<nbSamples;i++){
+ t = tMin + i * DT;
+ const SVector3 dp = q - position(t);
+ const double D = dp.norm();
+ if (D < DMIN) {
+ topt = t;
+ DMIN = D;
}
- if (i > nb) relax *= 0.5;
}
+
+ // printf("parameter %g as an initial guess (dist = %g)\n",topt,DMIN);
+
+ if (topt == tMin)
+ t = goldenSectionSearch (this, q, topt, topt + DT/2, topt + DT, 1.e-7);
+ else if (topt == tMax)
+ t = goldenSectionSearch (this, q, topt - DT, topt - DT/2 , topt, 1.e-7);
+ else
+ t = goldenSectionSearch (this, q, topt - DT, topt, topt + DT, 1.e-7);
+
+ const SVector3 dp = q - position(t);
+ const double D = dp.norm();
+
+ // printf("after golden section parameter %g (dist = %g)\n",t,D);
+
return point(t);
}
@@ -348,13 +393,13 @@ bool GEdge::XYZToU(const double X, const double Y, const double Z,
}
if(relax > 1.e-2) {
- Msg::Info("point %g %g %g on edge %d : Relaxation factor = %g",
- Q.x(), Q.y(), Q.z(), 0.75 * relax);
+ // Msg::Info("point %g %g %g on edge %d : Relaxation factor = %g",
+ // Q.x(), Q.y(), Q.z(), 0.75 * relax);
return XYZToU(Q.x(), Q.y(), Q.z(), u, 0.75 * relax);
}
- Msg::Error("Could not converge reparametrisation of point (%e,%e,%e) on edge %d",
- Q.x(), Q.y(), Q.z(), tag());
+ // Msg::Error("Could not converge reparametrisation of point (%e,%e,%e) on edge %d",
+ // Q.x(), Q.y(), Q.z(), tag());
return false;
}
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index ce18e3bb22..06b6889550 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -228,6 +228,7 @@ GEdge *OCCFactory::addCircleArc(GModel *gm, const arcCreationMethod &method,
gp_Pnt aP3(end->x(), end->y(), end->z());
TopoDS_Edge occEdge;
+
OCCVertex *occv1 = dynamic_cast<OCCVertex*>(start);
OCCVertex *occv2 = dynamic_cast<OCCVertex*>(end);
@@ -241,7 +242,7 @@ GEdge *OCCFactory::addCircleArc(GModel *gm, const arcCreationMethod &method,
}
else if (method == GModelFactory::CENTER_START_END){
Standard_Real Radius = aP1.Distance(aP2);
- gce_MakeCirc MC(aP1,gce_MakePln(aP1, aP2, aP3).Value(), Radius);
+ gce_MakeCirc MC(aP2,gce_MakePln(aP1, aP2, aP3).Value(), Radius);
const gp_Circ& Circ = MC.Value();
Standard_Real Alpha1 = ElCLib::Parameter(Circ, aP1);
Standard_Real Alpha2 = ElCLib::Parameter(Circ, aP3);
diff --git a/Geo/GeomMeshMatcher.cpp b/Geo/GeomMeshMatcher.cpp
index 04b501d11d..2d524a234f 100644
--- a/Geo/GeomMeshMatcher.cpp
+++ b/Geo/GeomMeshMatcher.cpp
@@ -17,6 +17,8 @@
#include "GmshMessage.h"
#include "SOrientedBoundingBox.h"
#include "MElement.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
#include "MVertex.h"
#include "MLine.h"
#include "MPoint.h"
@@ -423,6 +425,173 @@ void GeomMeshMatcher::destroy()
delete GeomMeshMatcher::_gmm_instance;
}
+static GVertex *getGVertex (MVertex *v1, GModel *gm, const double TOL){
+ GVertex *best = 0;
+ double bestScore = TOL;
+ for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
+ {
+ GVertex *v2 = (*it)->getBeginVertex();
+ double score = sqrt((v1->x() - v2->x())*(v1->x() - v2->x()) +
+ (v1->y() - v2->y())*(v1->y() - v2->y()) +
+ (v1->z() - v2->z())*(v1->z() - v2->z()));
+ if (score < bestScore){
+ bestScore = score;
+ best = v2;
+ }
+ }
+ {
+ GVertex *v2 = (*it)->getEndVertex();
+ double score = sqrt((v1->x() - v2->x())*(v1->x() - v2->x()) +
+ (v1->y() - v2->y())*(v1->y() - v2->y()) +
+ (v1->z() - v2->z())*(v1->z() - v2->z()));
+ if (score < bestScore){
+ bestScore = score;
+ best = v2;
+ }
+ }
+ }
+ // if (best)printf("getting point %g %g on vertices best score is %12.5E\n",v1->x(),v1->y(),bestScore);
+ return best;
+}
+
+static GPoint getGEdge (MVertex *v1, GModel *gm, const double TOL){
+ GPoint gpBest;
+ double bestScore = TOL;
+
+
+ for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
+ GEdge *e = *it;
+ double pp;
+ GPoint gp = e->closestPoint(SPoint3(v1->x(),v1->y(),v1->z()), pp);
+ double score = sqrt((v1->x() - gp.x())*(v1->x() - gp.x()) +
+ (v1->y() - gp.y())*(v1->y() - gp.y()) +
+ (v1->z() - gp.z())*(v1->z() - gp.z()));
+ if (score < bestScore){
+ bestScore = score;
+ gpBest = gp;
+ }
+ }
+
+ // printf("getting point %g %g (%g %g) on edges best score is %12.5E\n",v1->x(),v1->y(),gpBest.x(),gpBest.y(),bestScore);
+ return gpBest;
+}
+
+static GPoint getGFace (MVertex *v1, GModel *gm, const double TOL){
+ GPoint gpBest;
+ double bestScore = TOL;
+ for (GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){
+ GFace *gf = *it;
+ SPoint2 pp;
+ double guess[2] = {0,0};
+ GPoint gp = gf->closestPoint(SPoint3(v1->x(),v1->y(),v1->z()), guess);
+ double score = sqrt((v1->x() - gp.x())*(v1->x() - gp.x()) +
+ (v1->y() - gp.y())*(v1->y() - gp.y()) +
+ (v1->z() - gp.z())*(v1->z() - gp.z()));
+ if (score < bestScore){
+ bestScore = score;
+ gpBest = gp;
+ }
+ }
+ // printf("getting point %g %g (%g %g) on faces best score is %12.5E\n",v1->x(),v1->y(),gpBest.x(),gpBest.y(),bestScore);
+ return gpBest;
+}
+
+
+int GeomMeshMatcher::forceTomatch(GModel *geom, GModel *mesh, const double TOL)
+{
+ // assume that the geometry is the right one
+
+ std::vector<GEntity*> entities;
+ mesh->getEntities(entities);
+ for(unsigned int i = 0; i < entities.size(); i++){
+ for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++){
+ MVertex *v = entities[i]->mesh_vertices[j];
+ GVertex *gv = getGVertex (v, geom, TOL);
+ bool found = 0;
+ if (gv){
+ printf("vertex %d matches GVertex %d\n",v->getNum(),gv->tag());
+ found=1;
+ MVertex *vvv = new MVertex (v->x(),v->y(),v->z(),gv,v->getNum());
+ gv->mesh_vertices.push_back(vvv);
+ gv->points.push_back(new MPoint(vvv,v->getNum()));
+
+ }
+ else if (v->onWhat()->dim() == 1){
+ GPoint gp = getGEdge (v, geom, 1.e22);
+ if(gp.g()){
+ GEntity *gg = (GEntity*)gp.g();
+ found=1;
+ // printf("vertex %d matches GEdge %d on position %g\n",v->getNum(),gg->tag(),gp.u());
+ gg->mesh_vertices.push_back(new MEdgeVertex (gp.x(),gp.y(),gp.z(),
+ gg,gp.u(),-1.,v->getNum()));
+ }
+ }
+ if (!found && v->onWhat()->dim() <= 2){
+ GPoint gp = getGFace (v, geom, TOL);
+ if(gp.g()){
+ GEntity *gg = (GEntity*)gp.g();
+ found = 1;
+ // printf("vertex %d matches GFace %d\n",v->getNum(),gg->tag());
+ gg->mesh_vertices.push_back(new MFaceVertex (gp.x(),gp.y(),gp.z(),
+ gg,gp.u(),gp.v(),v->getNum()));
+ }
+ }
+ if (!found) Msg::Error("vertex %d classified on %d %d not matched",v->getNum(),v->onWhat()->dim(),v->onWhat()->tag());
+ }
+ }
+ // printf("creating edges\n");
+ for (GModel::eiter it = mesh->firstEdge(); it != mesh->lastEdge(); ++it){
+ // printf("edge %d\n",(*it)->tag());
+ for (int i=0;i<(*it)->lines.size();i++){
+ // printf("medge %d %d\n",(*it)->lines[i]->getVertex(0)->getNum(),(*it)->lines[i]->getVertex(1)->getNum());
+ MVertex *v1 = geom->getMeshVertexByTag((*it)->lines[i]->getVertex(0)->getNum());
+ MVertex *v2 = geom->getMeshVertexByTag((*it)->lines[i]->getVertex(1)->getNum());
+ if (v1 && v2){
+ GEdge *ge= 0;
+ if (v1->onWhat()->dim() == 1)ge = (GEdge*)v1->onWhat();
+ if (v2->onWhat()->dim() == 1)ge = (GEdge*)v2->onWhat();
+ if (ge){
+ double u1,u2;
+ reparamMeshVertexOnEdge(v1, ge, u1);
+ reparamMeshVertexOnEdge(v2, ge, u2);
+ if (u1< u2)ge->lines.push_back(new MLine(v1,v2));
+ else ge->lines.push_back(new MLine(v2,v1));
+ }
+ else printf("argh !\n");
+ }
+ else{
+ if (!v1)printf("Vertex %d has not been found\n", (*it)->lines[i]->getVertex(0)->getNum());
+ if (!v2)printf("Vertex %d has not been found\n", (*it)->lines[i]->getVertex(1)->getNum());
+ }
+ }
+ }
+ printf("creating faces\n");
+ for (GModel::fiter it = mesh->firstFace(); it != mesh->lastFace(); ++it){
+ for (int i=0;i<(*it)->triangles.size();i++){
+ MVertex *v1 = geom->getMeshVertexByTag((*it)->triangles[i]->getVertex(0)->getNum());
+ MVertex *v2 = geom->getMeshVertexByTag((*it)->triangles[i]->getVertex(1)->getNum());
+ MVertex *v3 = geom->getMeshVertexByTag((*it)->triangles[i]->getVertex(2)->getNum());
+ if (v1->onWhat()->dim() == 2)((GFace*)v1->onWhat())->triangles.push_back(new MTriangle(v1,v2,v3));
+ else if (v2->onWhat()->dim() == 2)((GFace*)v2->onWhat())->triangles.push_back(new MTriangle(v1,v2,v3));
+ else if (v3->onWhat()->dim() == 2)((GFace*)v3->onWhat())->triangles.push_back(new MTriangle(v1,v2,v3));
+ }
+ for (int i=0;i<(*it)->quadrangles.size();i++){
+ MVertex *v1 = geom->getMeshVertexByTag((*it)->quadrangles[i]->getVertex(0)->getNum());
+ MVertex *v2 = geom->getMeshVertexByTag((*it)->quadrangles[i]->getVertex(1)->getNum());
+ MVertex *v3 = geom->getMeshVertexByTag((*it)->quadrangles[i]->getVertex(2)->getNum());
+ MVertex *v4 = geom->getMeshVertexByTag((*it)->quadrangles[i]->getVertex(3)->getNum());
+ // printf("quad %p %p %p %p\n",v1,v2,v3,v4);
+ if (v1->onWhat()->dim() == 2)((GFace*)v1->onWhat())->quadrangles.push_back(new MQuadrangle(v1,v2,v3,v4));
+ else if (v2->onWhat()->dim() == 2)((GFace*)v2->onWhat())->quadrangles.push_back(new MQuadrangle(v1,v2,v3,v4));
+ else if (v3->onWhat()->dim() == 2)((GFace*)v3->onWhat())->quadrangles.push_back(new MQuadrangle(v1,v2,v3,v4));
+ else if (v4->onWhat()->dim() == 2)((GFace*)v4->onWhat())->quadrangles.push_back(new MQuadrangle(v1,v2,v3,v4));
+ }
+ }
+ geom->writeMSH("hopla.msh",2.2,false,false,true);
+}
+
+
+
int GeomMeshMatcher::match(GModel *geom, GModel *mesh)
{
mesh->createTopologyFromMesh();
diff --git a/Geo/GeomMeshMatcher.h b/Geo/GeomMeshMatcher.h
index 7ec6f51002..a05eff3dd7 100644
--- a/Geo/GeomMeshMatcher.h
+++ b/Geo/GeomMeshMatcher.h
@@ -35,6 +35,8 @@ class GeomMeshMatcher {
static GeomMeshMatcher *instance();
static void destroy();
int match(GModel* geom, GModel* mesh);
+ int forceTomatch(GModel *geom, GModel *mesh, const double TOL);
+
};
#endif
diff --git a/Geo/MElement.h b/Geo/MElement.h
index b4836a4dc1..534f22c1a4 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -246,6 +246,16 @@ class MElement
// parametric coordinates
double getJacobian(const fullMatrix<double> &gsf, double jac[3][3]);
double getJacobian(double u, double v, double w, double jac[3][3]);
+ inline double getJacobian(double u, double v, double w, fullMatrix<double> &j){
+ double JAC[3][3];
+ const double detJ = getJacobian (u,v,w,JAC);
+ for (int i=0;i<3;i++){
+ j(i,0) = JAC[i][0];
+ j(i,1) = JAC[i][1];
+ j(i,2) = JAC[i][2];
+ }
+ return detJ;
+ }
double getPrimaryJacobian(double u, double v, double w, double jac[3][3]);
double getJacobianDeterminant(double u, double v, double w)
{
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 07ef9fcd35..22aeb98b05 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -23,6 +23,7 @@ set(SRC
BDS.cpp
HighOrder.cpp
highOrderSmoother.cpp
+ highOrderTools.cpp
meshPartition.cpp
meshPartitionOptions.cpp
meshRefine.cpp
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index b34adb5348..e03cc759a8 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -9,6 +9,7 @@
#include "HighOrder.h"
#include "highOrderSmoother.h"
+#include "highOrderTools.h"
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
@@ -256,7 +257,7 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
if(reparamOK){
double relax = 1.;
while (1){
- if(computeEquidistantParameters(ge, u0, u1, nPts + 2, US, relax))
+ if(computeEquidistantParameters(ge, std::min(u0,u1), std::max(u0,u1), nPts + 2, US, relax))
break;
relax /= 2.0;
if(relax < 1.e-2)
@@ -264,8 +265,8 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
}
if(relax < 1.e-2)
Msg::Warning
- ("Failed to compute equidistant parameters (relax = %g) for edge %d-%d",
- relax, v0->getNum(), v1->getNum());
+ ("Failed to compute equidistant parameters (relax = %g, value = %g) for edge %d-%d parametrized with %g %g on GEdge %d",
+ relax, US[1], v0->getNum(), v1->getNum(),u0,u1,ge->tag());
}
}
for(int j = 0; j < nPts; j++){
@@ -273,7 +274,7 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
double uc = (1. - t) * u0 + t * u1; // can be wrong, that's ok
MVertex *v;
- if(linear || !reparamOK || uc < u0 || uc > u1){
+ if(linear || !reparamOK || uc < std::min(u0,u1) || uc > std::max(u0,u1)){
Msg::Warning("We don't have a valid parameter on curve %d-%d",
v0->getNum(), v1->getNum());
// we don't have a (valid) parameter on the curve
@@ -281,8 +282,8 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
v = new MVertex(pc.x(), pc.y(), pc.z(), ge);
}
else {
- GPoint pc = ge->point(US[j + 1]);
- v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge, US[j + 1]);
+ GPoint pc = ge->point(US[u0<u1? j + 1 : nPts - 1 - (j + 1)]);
+ v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge, US[u0<u1? j + 1 : nPts - 1 - (j + 1)]);
if (displ2D || displ3D){
SPoint3 pc2 = edge.interpolate(t);
@@ -1212,10 +1213,10 @@ void checkHighOrderTriangles(const char* cc, GModel *m,
}
}
if (minJGlob > 0)
- Msg::Info("%s : Worst Face Smoothness %g Gamma %g NbFair = %d",
+ Msg::Info("%s : Worst Face Distorsion Mapping %g Gamma %g Nb elem. (0<d<0.2) = %d",
cc, minJGlob, minGGlob,nbfair );
else
- Msg::Warning("%s : Worst Face Smoothness %g (%d negative jacobians) "
+ Msg::Warning("%s : Worst Face Distorsion Mapping %g (%d negative jacobians) "
"Worst Gamma %g Avg Smoothness %g", cc, minJGlob, bad.size(),
minGGlob, avg / (count ? count : 1));
}
@@ -1252,7 +1253,7 @@ static void checkHighOrderTetrahedron(const char* cc, GModel *m,
extern double mesh_functional_distorsion(MElement *t, double u, double v);
-static void printJacobians(GModel *m, const char *nm)
+void printJacobians(GModel *m, const char *nm)
{
const int n = 100;
double D[n][n], X[n][n], Y[n][n], Z[n][n];
@@ -1380,5 +1381,10 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// printJacobians(m, "smoothness.pos");
double t2 = Cpu();
+
+ highOrderTools hot (m);
+ hot.ensureMinimumDistorsion(0.1);
+ checkHighOrderTriangles("final mesh", m, bad, worst);
+
Msg::StatusBar(2, true, "Done meshing order %d (%g s)", order, t2 - t1);
}
diff --git a/Mesh/highOrderSmoother.cpp b/Mesh/highOrderSmoother.cpp
index 6992b2f64b..42f21fcbf2 100644
--- a/Mesh/highOrderSmoother.cpp
+++ b/Mesh/highOrderSmoother.cpp
@@ -483,7 +483,7 @@ void highOrderSmoother::smooth_metric(std::vector<MElement*> & all, GFace *gf)
std::vector<MElement*> layer, v;
double minD;
- getDistordedElements(all, 0.5, v, minD);
+ getDistordedElements(all, 0.99, v, minD);
const int nbLayers = 3; //2, originally :)
for (int i = 0; i < nbLayers; i++){
@@ -591,11 +591,6 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
if (myAssembler.sizeOfR()){
- // initialize _materialLaw to 1 everywhere
- for (unsigned int i = 0; i < v.size(); i++){
- MElement *e = v[i];
- // _materialLaw[e] = 1.0;
- }
// while convergence -------------------------------------------------------
for (unsigned int i = 0; i < v.size(); i++){
MElement *e = v[i];
@@ -612,8 +607,6 @@ double highOrderSmoother::smooth_metric_(std::vector<MElement*> & v,
fullMatrix<double> J23K33(n2, n3);
K33.setAll(0.0);
SElement se(e);
- // set kappa to the actual kappa of the material law
- // El.setCompressibility( compressibilityFunction ( _materialLaw[e] ) );
El.elementMatrix(&se, K33);
computeMetricVector(gf, e, J32, J23, D3);
J23K33.gemm(J23, K33, 1, 0);
@@ -861,7 +854,7 @@ void highOrderSmoother::smooth(std::vector<MElement*> &all)
std::vector<MElement*> layer, v;
double minD;
- getDistordedElements(all, 0.5, v, minD);
+ getDistordedElements(all, 0.85, v, minD);
Msg::Info("%d elements / %d distorted min Disto = %g\n",
all.size(), v.size(), minD);
@@ -1756,5 +1749,71 @@ void highOrderSmoother::smooth_pNpoint(GFace *gf)
findOptimalLocationsPN(gf,this);
}
+//-------------------------------------------------
+// Assume a GModel that is meshed with high order nodes that
+// are on their final position (at least they all lie on
+// surfaces and curves and the smoother may move them)
+//-------------------------------------------------
+
+highOrderSmoother::highOrderSmoother (GModel *gm)
+ : _tag(111)
+{
+ _clean();
+ // compute straigh sided positions that are actually X,Y,Z positions
+ // that are NOT always on curves and surfaces
+
+ // compute stright sided positions of vertices that are classified on model edges
+ for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
+ for (int i=0;i<(*it)->lines.size();i++){
+ MLine *l = (*it)->lines[i];
+ int N = l->getNumVertices()-2;
+ SVector3 p0((*it)->lines[i]->getVertex(0)->x(),
+ (*it)->lines[i]->getVertex(0)->y(),
+ (*it)->lines[i]->getVertex(0)->z());
+ SVector3 p1((*it)->lines[i]->getVertex(1)->x(),
+ (*it)->lines[i]->getVertex(1)->y(),
+ (*it)->lines[i]->getVertex(1)->z());
+ for (int j=1;j<=N;j++){
+ const double xi = (double)j/(N+1);
+ const SVector3 straightSidedPoint = p0 *(1.-xi) + p1*xi;
+ _straightSidedLocation [(*it)->lines[i]->getVertex(j)] = straightSidedPoint;
+ }
+ }
+ }
+
+ // compute stright sided positions of vertices that are classified on model faces
+ for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){
+ for (int i=0;i<(*it)->triangles.size();i++){
+ MTriangle *t = (*it)->triangles[i];
+ MFace face = t->getFace(0);
+ const polynomialBasis* fs = t->getFunctionSpace();
+ for (int j=0;j<t->getNumVertices();j++){
+ if (t->getVertex(j)->onWhat() == *it){
+ const double t1 = fs->points(j, 0);
+ const double t2 = fs->points(j, 1);
+ SPoint3 pc = face.interpolate(t1, t2);
+ _straightSidedLocation [t->getVertex(j)] =
+ SVector3(pc.x(),pc.y(),pc.z());
+ }
+ }
+ }
+ for (int i=0;i<(*it)->quadrangles.size();i++){
+ MQuadrangle *q = (*it)->quadrangles[i];
+ MFace face = q->getFace(0);
+ const polynomialBasis* fs = q->getFunctionSpace();
+ for (int j=0;j<q->getNumVertices();j++){
+ if (q->getVertex(j)->onWhat() == *it){
+ const double t1 = fs->points(j, 0);
+ const double t2 = fs->points(j, 1);
+ SPoint3 pc = face.interpolate(t1, t2);
+ _straightSidedLocation [q->getVertex(j)] =
+ SVector3(pc.x(),pc.y(),pc.z());
+ }
+ }
+ }
+ }
+}
+
+
#endif
diff --git a/Mesh/highOrderSmoother.h b/Mesh/highOrderSmoother.h
index ef1303ef89..3c5439c822 100644
--- a/Mesh/highOrderSmoother.h
+++ b/Mesh/highOrderSmoother.h
@@ -26,7 +26,6 @@ class GRegion;
class highOrderSmoother
{
const int _tag;
- std::map<MElement*,std::vector<double> > _materialLaw;
std::map<MVertex*,SVector3> _straightSidedLocation;
std::map<MVertex*,SVector3> _targetLocation;
int _dim;
diff --git a/Mesh/highOrderTools.cpp b/Mesh/highOrderTools.cpp
new file mode 100644
index 0000000000..5c8df2bf57
--- /dev/null
+++ b/Mesh/highOrderTools.cpp
@@ -0,0 +1,668 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Koen Hillewaert
+//
+
+#include "GmshConfig.h"
+
+#if defined(HAVE_SOLVER)
+
+#include "MLine.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+#include "MTetrahedron.h"
+#include "MHexahedron.h"
+#include "MPrism.h"
+#include "MPyramid.h"
+#include "MPoint.h"
+#include "HighOrder.h"
+#include "meshGFaceOptimize.h"
+#include "highOrderTools.h"
+#include "GFace.h"
+#include "GRegion.h"
+#include "GeomMeshMatcher.h"
+#include "Context.h"
+#include "Numeric.h"
+#include "dofManager.h"
+#include "elasticityTerm.h"
+#include "linearSystemCSR.h"
+#include "linearSystemFull.h"
+#include "linearSystemPETSc.h"
+
+#define SQU(a) ((a)*(a))
+
+void highOrderTools::moveToStraightSidedLocation(MElement *e) const
+{
+ for(int i = 0; i < e->getNumVertices(); i++){
+ MVertex *v = e->getVertex(i);
+ std::map<MVertex*,SVector3>::const_iterator it = _straightSidedLocation.find(v);
+ if (it != _straightSidedLocation.end()){
+ // printf("move %d : %g %g -> %g %g\n",v->getNum(),v->x() ,v->y() , it->second.x() , it->second.y() );
+ v->x() = it->second.x();
+ v->y() = it->second.y();
+ v->z() = it->second.z();
+ }
+ }
+}
+
+void highOrderTools::ensureMinimumDistorsion(MElement *e, double threshold)
+{
+ double dist = e->distoShapeMeasure();
+ if (dist > threshold) return;
+ double a1 = 0., a2 = 1.0, x[1000][3], X[1000][3] ;
+ for(int i = 0; i < e->getNumVertices(); i++){
+ MVertex *v = e->getVertex(i);
+ x[i][0] = v->x();
+ x[i][1] = v->y();
+ x[i][2] = v->z();
+ std::map<MVertex*,SVector3>::const_iterator it = _straightSidedLocation.find(v);
+ if (it != _straightSidedLocation.end()){
+ X[i][0] = it->second.x();
+ X[i][1] = it->second.y();
+ X[i][2] = it->second.z();
+ }
+ else {
+ X[i][0] = v->x();
+ X[i][1] = v->y();
+ X[i][2] = v->z();
+ }
+ }
+ // a == 0 -> straight
+ // a == 1 -> curved
+ while(1){
+ double a = 0.5*(a1+a2);
+ for(int i = 0; i < e->getNumVertices(); i++){
+ MVertex *v = e->getVertex(i);
+ v->x() = a * x[i][0] + (1.-a) * X[i][0];
+ v->y() = a * x[i][1] + (1.-a) * X[i][1];
+ v->z() = a * x[i][2] + (1.-a) * X[i][2];
+ }
+ double dist = e->distoShapeMeasure();
+ // printf("a = %g dist = %g\n",a,dist);
+ // getchar();
+ if (dist > 0 && fabs(dist-threshold) < .05)break;
+ if (dist < threshold)a2 = a;
+ else a1 = a;
+ }
+ // printf("element done\n");
+}
+
+
+static void getDistordedElements(const std::vector<MElement*> &v,
+ const double & threshold,
+ std::vector<MElement*> &d,
+ double &minD)
+{
+ d.clear();
+ minD = 1;
+ for (unsigned int i = 0; i < v.size(); i++){
+ const double disto = v[i]->distoShapeMeasure();
+ if (disto < threshold)
+ d.push_back(v[i]);
+ minD = std::min(minD, disto);
+ }
+}
+
+static void addOneLayer(const std::vector<MElement*> &v,
+ std::vector<MElement*> &d,
+ std::vector<MElement*> &layer)
+{
+ std::set<MVertex*> all;
+ for (unsigned int i = 0; i < d.size(); i++){
+ MElement *e = d[i];
+ int n = e->getNumPrimaryVertices();
+ for (int j = 0; j < n; j++){
+ all.insert(e->getVertex(j));
+ }
+ }
+ layer.clear();
+ std::sort(d.begin(), d.end());
+
+ for (unsigned int i = 0; i < v.size(); i++){
+ MElement *e = v[i];
+ bool found = std::binary_search(d.begin(), d.end(), e);
+ // element is not yet there
+ if (!found){
+ int n = e->getNumPrimaryVertices();
+ for (int j = 0; j < n; j++){
+ MVertex *vert = e->getVertex(j);
+ if (all.find(vert) != all.end()){
+ layer.push_back(e);
+ j = n;
+ }
+ }
+ }
+ }
+}
+
+
+double highOrderTools::applySmoothingTo(GFace *gf, double tres, bool mixed)
+{
+ if (!gf){
+ Msg::Error("Cannot smooth that face");
+ return -1;
+ }
+ std::vector<MElement*> v;
+ v.insert(v.begin(), gf->triangles.begin(),gf->triangles.end());
+ v.insert(v.begin(), gf->quadrangles.begin(),gf->quadrangles.end());
+ Msg::Info("Smoothing high order mesh : model face %d (%d elements)", gf->tag(),
+ v.size());
+ // applySmoothingTo(v,gf);
+ return applySmoothingTo(v,tres, mixed);
+}
+
+
+void highOrderTools::ensureMinimumDistorsion (double threshold){
+ std::vector<MElement*> v;
+ if (_dim == 2){
+ for (GModel::fiter fit = _gm->firstFace(); fit != _gm->lastFace(); ++fit) {
+ v.insert(v.begin(), (*fit)->triangles.begin(),(*fit)->triangles.end());
+ v.insert(v.begin(), (*fit)->quadrangles.begin(),(*fit)->quadrangles.end());
+ }
+ }
+ else if (_dim == 3){
+ for (GModel::riter rit = _gm->firstRegion(); rit != _gm->lastRegion(); ++rit) {
+ v.insert(v.begin(), (*rit)->hexahedra.begin(),(*rit)->hexahedra.end());
+ v.insert(v.begin(), (*rit)->tetrahedra.begin(),(*rit)->tetrahedra.end());
+ v.insert(v.begin(), (*rit)->prisms.begin(),(*rit)->prisms.end());
+ }
+ }
+ ensureMinimumDistorsion(v,threshold);
+}
+
+
+/*
+void highOrderTools::applySmoothingTo(GRegion *gr)
+{
+ std::vector<MElement*> v;
+ v.insert(v.begin(), gr->tetrahedra.begin(),gr->tetrahedra.end());
+ v.insert(v.begin(), gr->hexahedra.begin(),gr->hexahedra.end());
+ v.insert(v.begin(), gr->prisms.begin(),gr->prisms.end());
+ Msg::Info("Smoothing high order mesh : model region %d (%d elements)", gr->tag(),
+ v.size());
+ applySmoothingTo(v,gr);
+}
+*/
+
+void highOrderTools::computeMetricInfo(GFace *gf,
+ MElement *e,
+ fullMatrix<double> &J,
+ fullMatrix<double> &JT,
+ fullVector<double> &D)
+{
+ int nbNodes = e->getNumVertices();
+ for (int j = 0; j < nbNodes; j++){
+ SPoint2 param;
+ reparamMeshVertexOnFace(e->getVertex(j), gf, param);
+ Pair<SVector3,SVector3> der = gf->firstDer(param);
+ int XJ = j;
+ int YJ = j + nbNodes;
+ int ZJ = j + 2 * nbNodes;
+ int UJ = j;
+ int VJ = j + nbNodes;
+ J(XJ,UJ) = der.first().x();
+ J(YJ,UJ) = der.first().y();
+ J(ZJ,UJ) = der.first().z();
+ J(XJ,VJ) = der.second().x();
+ J(YJ,VJ) = der.second().y();
+ J(ZJ,VJ) = der.second().z();
+
+ JT(UJ,XJ) = der.first().x();
+ JT(UJ,YJ) = der.first().y();
+ JT(UJ,ZJ) = der.first().z();
+ JT(VJ,XJ) = der.second().x();
+ JT(VJ,YJ) = der.second().y();
+ JT(VJ,ZJ) = der.second().z();
+
+ GPoint gp = gf->point(param);
+ SVector3 ss = getSSL(e->getVertex(j));
+ D(XJ) = gp.x() - ss.x();
+ D(YJ) = gp.y() - ss.y();
+ D(ZJ) = gp.z() - ss.z();
+ }
+}
+
+void highOrderTools::applySmoothingTo(std::vector<MElement*> & all, GFace *gf)
+{
+#ifdef HAVE_TAUCS
+ linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
+#else
+ linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
+#endif
+
+ dofManager<double> myAssembler(lsys);
+ elasticityTerm El(0, 1.0, .33, _tag);
+ std::vector<MElement*> layer, v;
+ double minD;
+ getDistordedElements(all, 0.5, v, minD);
+ const int nbLayers = 3;
+ for (int i = 0; i < nbLayers; i++){
+ addOneLayer(all, v, layer);
+ v.insert(v.end(), layer.begin(), layer.end());
+ }
+ Msg::Debug("%d elements after adding %d layers", (int)v.size(), nbLayers);
+
+ if (!v.size()) return;
+ addOneLayer(all, v, layer);
+ std::set<MVertex*>::iterator it;
+ std::map<MVertex*,SVector3>::iterator its;
+ std::map<MVertex*,SVector3>::iterator itpresent;
+ std::map<MVertex*,SVector3>::iterator ittarget;
+ std::set<MVertex*> verticesToMove;
+
+ // on the last layer, fix displacement to 0
+ for (unsigned int i = 0; i < layer.size(); i++){
+ for (int j = 0; j < layer[i]->getNumVertices(); j++){
+ MVertex *vert = layer[i]->getVertex(j);
+ myAssembler.fixVertex(vert, 0, _tag, 0);
+ myAssembler.fixVertex(vert, 1, _tag, 0);
+ }
+ }
+
+ // fix all vertices that cannot move
+ for (unsigned int i = 0; i < v.size(); i++){
+ for (int j = 0; j < v[i]->getNumVertices(); j++){
+ MVertex *vert = v[i]->getVertex(j);
+ if (vert->onWhat()->dim() < 2){
+ myAssembler.fixVertex(vert, 0, _tag, 0);
+ myAssembler.fixVertex(vert, 1, _tag, 0);
+ }
+ }
+ }
+
+ // number the other DOFs
+ for (unsigned int i = 0; i < v.size(); i++){
+ for (int j = 0; j < v[i]->getNumVertices(); j++){
+ MVertex *vert = v[i]->getVertex(j);
+ myAssembler.numberVertex(vert, 0, _tag);
+ myAssembler.numberVertex(vert, 1, _tag);
+ verticesToMove.insert(vert);
+ }
+ }
+
+ double dx0 = smooth_metric_(v, gf, myAssembler, verticesToMove, El);
+ double dx = dx0;
+ Msg::Debug(" dx0 = %12.5E", dx0);
+ int iter = 0;
+ while(1){
+ double dx2 = smooth_metric_(v, gf, myAssembler, verticesToMove, El);
+ Msg::Debug(" dx2 = %12.5E", dx2);
+ if (fabs(dx2 - dx) < 1.e-4 * dx0)break;
+ if (iter++ > 2)break;
+ dx = dx2;
+ }
+
+ for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
+ SPoint2 param;
+ reparamMeshVertexOnFace(*it, gf, param);
+ GPoint gp = gf->point(param);
+ if ((*it)->onWhat()->dim() == 2){
+ (*it)->x() = gp.x();
+ (*it)->y() = gp.y();
+ (*it)->z() = gp.z();
+ _targetLocation[*it] = SVector3(gp.x(), gp.y(), gp.z());
+ }
+ else{
+ SVector3 p = _targetLocation[(*it)];
+ (*it)->x() = p.x();
+ (*it)->y() = p.y();
+ (*it)->z() = p.z();
+ }
+ }
+ delete lsys;
+}
+
+double highOrderTools::smooth_metric_(std::vector<MElement*> & v,
+ GFace *gf,
+ dofManager<double> &myAssembler,
+ std::set<MVertex*> &verticesToMove,
+ elasticityTerm &El)
+{
+ std::set<MVertex*>::iterator it;
+ double dx = 0.0;
+
+ // printf("size %d\n",myAssembler.sizeOfR());
+
+ if (myAssembler.sizeOfR()){
+
+ // while convergence -------------------------------------------------------
+ for (unsigned int i = 0; i < v.size(); i++){
+ MElement *e = v[i];
+ int nbNodes = e->getNumVertices();
+ const int n2 = 2 * nbNodes;
+ const int n3 = 3 * nbNodes;
+
+ fullMatrix<double> K33(n3, n3);
+ fullMatrix<double> K22(n2, n2);
+ fullMatrix<double> J32(n3, n2);
+ fullMatrix<double> J23(n2, n3);
+ fullVector<double> D3(n3);
+ fullVector<double> R2(n2);
+ fullMatrix<double> J23K33(n2, n3);
+ K33.setAll(0.0);
+ SElement se(e);
+ El.elementMatrix(&se, K33);
+ computeMetricInfo(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++){
+ Dof RDOF = El.getLocalDofR(&se, j);
+ myAssembler.assemble(RDOF, -R2(j));
+ for (int k = 0; k < n2; k++){
+ Dof CDOF = El.getLocalDofC(&se, k);
+ myAssembler.assemble(RDOF, CDOF, K22(j, k));
+ }
+ }
+ }
+ myAssembler.systemSolve();
+ // for all element, compute detJ at integration points --> material law
+ // end while convergence -------------------------------------------------------
+
+ for (it = verticesToMove.begin(); it != verticesToMove.end(); ++it){
+ if ((*it)->onWhat()->dim() == 2){
+ SPoint2 param;
+ reparamMeshVertexOnFace((*it), gf, param);
+ SPoint2 dparam;
+ myAssembler.getDofValue((*it), 0, _tag, dparam[0]);
+ myAssembler.getDofValue((*it), 1, _tag, dparam[1]);
+ SPoint2 newp = param+dparam;
+ dx += newp.x() * newp.x() + newp.y() * newp.y();
+ (*it)->setParameter(0, newp.x());
+ (*it)->setParameter(1, newp.y());
+ }
+ }
+ myAssembler.systemClear();
+ }
+
+ return dx;
+}
+
+
+highOrderTools::highOrderTools (GModel *gm) : _gm(gm), _dim(2), _tag(111)
+{
+ computeStraightSidedPositions();
+}
+
+highOrderTools::highOrderTools (GModel *gm, GModel *mesh, int order) : _gm(gm), _dim(2), _tag(111)
+{
+ GeomMeshMatcher::instance()->forceTomatch(gm,mesh,1.e-5);
+ GeomMeshMatcher::instance()->destroy();
+ SetOrderN(gm, order, false, false);
+ computeStraightSidedPositions();
+}
+
+
+void highOrderTools::computeStraightSidedPositions ()
+{
+ _clean();
+ // compute straigh sided positions that are actually X,Y,Z positions
+ // that are NOT always on curves and surfaces
+
+ // points classified on model vertices shall not move !
+ for(GModel::viter it = _gm->firstVertex(); it != _gm->lastVertex(); ++it){
+ if ((*it)->points.size()){
+ MPoint *p = (*it)->points[0];
+ MVertex *v = p->getVertex(0);
+ _straightSidedLocation [v] = SVector3((*it)->x(),(*it)->y(),(*it)->z());
+ _targetLocation [v] = SVector3((*it)->x(),(*it)->y(),(*it)->z());
+ }
+ }
+
+ // printf("coucou2\n");
+ // compute stright sided positions of vertices that are classified on model edges
+ for(GModel::eiter it = _gm->firstEdge(); it != _gm->lastEdge(); ++it){
+ for (int i=0;i<(*it)->lines.size();i++){
+ MLine *l = (*it)->lines[i];
+ int N = l->getNumVertices()-2;
+ SVector3 p0((*it)->lines[i]->getVertex(0)->x(),
+ (*it)->lines[i]->getVertex(0)->y(),
+ (*it)->lines[i]->getVertex(0)->z());
+ SVector3 p1((*it)->lines[i]->getVertex(1)->x(),
+ (*it)->lines[i]->getVertex(1)->y(),
+ (*it)->lines[i]->getVertex(1)->z());
+ for (int j=1;j<=N;j++){
+ const double xi = (double)j/(N+1);
+ // printf("xi = %g\n",xi);
+ const SVector3 straightSidedPoint = p0 *(1.-xi) + p1*xi;
+ MVertex *v = (*it)->lines[i]->getVertex(j+1);
+ if (_straightSidedLocation.find(v) == _straightSidedLocation.end()){
+ _straightSidedLocation [v] = straightSidedPoint;
+ _targetLocation[v] = SVector3(v->x(),v->y(),v->z());
+ }
+ }
+ }
+ }
+
+ // printf("coucou3\n");
+ // compute stright sided positions of vertices that are classified on model faces
+ for(GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it){
+ for (int i=0;i<(*it)->mesh_vertices.size();i++){
+ MVertex *v = (*it)->mesh_vertices[i];
+ _targetLocation[v] = SVector3(v->x(),v->y(),v->z());
+ }
+
+ for (int i=0;i<(*it)->triangles.size();i++){
+ MTriangle *t = (*it)->triangles[i];
+ MFace face = t->getFace(0);
+ const polynomialBasis* fs = t->getFunctionSpace();
+ for (int j=0;j<t->getNumVertices();j++){
+ if (t->getVertex(j)->onWhat() == *it){
+ const double t1 = fs->points(j, 0);
+ const double t2 = fs->points(j, 1);
+ SPoint3 pc = face.interpolate(t1, t2);
+ _straightSidedLocation [t->getVertex(j)] =
+ SVector3(pc.x(),pc.y(),pc.z());
+ }
+ }
+ }
+ for (int i=0;i<(*it)->quadrangles.size();i++){
+ // printf("coucou quad %d\n",i);
+ MQuadrangle *q = (*it)->quadrangles[i];
+ MFace face = q->getFace(0);
+ const polynomialBasis* fs = q->getFunctionSpace();
+ for (int j=0;j<q->getNumVertices();j++){
+ if (q->getVertex(j)->onWhat() == *it){
+ const double t1 = fs->points(j, 0);
+ const double t2 = fs->points(j, 1);
+ SPoint3 pc = face.interpolate(t1, t2);
+ _straightSidedLocation [q->getVertex(j)] =
+ SVector3(pc.x(),pc.y(),pc.z());
+ }
+ }
+ }
+ }
+
+ Msg::Info("highOrderTools has been set up : %d nodes are considered",_straightSidedLocation.size());
+}
+
+
+// apply a displacement that does not create elements that are
+// distorted over a value "thres"
+double highOrderTools::apply_incremental_displacement (double max_incr,
+ std::vector<MElement*> & v,
+ bool mixed,
+ double thres,
+ char *meshName,
+ std::vector<MElement*> & disto)
+{
+#ifdef HAVE_PETSC
+ // assume that the mesh is OK, yet already curved
+ linearSystemPETSc<double> *lsys = new linearSystemPETSc<double>;
+ lsys->setParameter("petscOptions","-pc_type ilu");
+
+ dofManager<double> myAssembler(lsys);
+ elasticityMixedTerm El_mixed (0, 1.0, .333, _tag);
+ elasticityTerm El (0, 1.0, .333, _tag);
+
+ std::set<MVertex*> _vertices;
+
+ //+++++++++ Boundary Conditions & Numbering +++++++++++++++++++++++++++++++
+ // fix all dof that correspond to vertices on the boundary
+ // the value is equal
+ for (unsigned int i = 0; i < v.size(); i++){
+ for (int j = 0; j < v[i]->getNumVertices(); j++){
+ MVertex *vert = v[i]->getVertex(j);
+ _vertices.insert(vert);
+ }
+ }
+
+ //+++++++++ Fix d tr(eps) = 0 +++++++++++++++++++++++++++++++
+ if (mixed){
+ for (unsigned int i = 0; i < disto.size(); i++){
+ for (int j = 0; j < disto[i]->getNumVertices(); j++){
+ MVertex *vert = disto[i]->getVertex(j);
+ myAssembler.fixVertex(vert, 4, _tag, 0.0);
+ }
+ }
+ }
+
+ for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end(); ++it){
+ MVertex *vert = *it;
+ std::map<MVertex*,SVector3>::iterator itt = _targetLocation.find(vert);
+ // impose displacement @ boundary
+ // if (!(vert->onWhat()->geomType() == GEntity::Line
+ // && vert->onWhat()->tag() == 2)){
+ if (itt != _targetLocation.end() && vert->onWhat()->dim() < _dim){
+ // printf("fixing motion of vertex %d to %g %g\n",vert->getNum(),itt->second.x()-vert->x(),itt->second.y()-vert->y());
+ myAssembler.fixVertex(vert, 0, _tag, itt->second.x()-vert->x());
+ myAssembler.fixVertex(vert, 1, _tag, itt->second.y()-vert->y());
+ myAssembler.fixVertex(vert, 2, _tag, itt->second.z()-vert->z());
+ }
+ // ensure we do not touch any vertex that is on the boundary
+ else if (vert->onWhat()->dim() < _dim){
+ myAssembler.fixVertex(vert, 0, _tag, 0);
+ myAssembler.fixVertex(vert, 1, _tag, 0);
+ myAssembler.fixVertex(vert, 2, _tag, 0);
+ }
+ // }
+ if (_dim == 2)myAssembler.fixVertex(vert, 2, _tag, 0);
+ // number vertices
+ myAssembler.numberVertex(vert, 0, _tag);
+ myAssembler.numberVertex(vert, 1, _tag);
+ myAssembler.numberVertex(vert, 2, _tag);
+ if (mixed){
+ myAssembler.numberVertex(vert, 3, _tag);
+ myAssembler.numberVertex(vert, 4, _tag);
+ }
+ }
+
+ //+++++++++ Assembly & Solve ++++++++++++++++++++++++++++++++++++
+ if (myAssembler.sizeOfR()){
+ // assembly of the elasticity term on the
+ for (unsigned int i = 0; i < v.size(); i++){
+ SElement se(v[i]);
+ if (mixed)El_mixed.addToMatrix(myAssembler, &se);
+ else El.addToMatrix(myAssembler, &se);
+ }
+ // solve the system
+ lsys->systemSolve();
+ }
+
+ //+++++++++ Move vertices @ maximum ++++++++++++++++++++++++++++++++++++
+ FILE *fd = fopen ("d.msh","w");
+ fprintf(fd,"$MeshFormat\n2 0 8\n$EndMeshFormat\n$NodeData\n1\n\"tr(sigma)\"\n1\n0.0\n3\n1\n3\n%d\n",_vertices.size());
+ for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end(); ++it){
+ double ax, ay, az;
+ myAssembler.getDofValue(*it, 0, _tag, ax);
+ myAssembler.getDofValue(*it, 1, _tag, ay);
+ myAssembler.getDofValue(*it, 2, _tag, az);
+ (*it)->x() += max_incr*ax;
+ (*it)->y() += max_incr*ay;
+ (*it)->z() += max_incr*az;
+ fprintf(fd,"%d %g %g %g\n",(*it)->getIndex(), ax,ay,az);
+ }
+ fprintf(fd,"$EndNodeData\n");
+ fclose(fd);
+
+ //+------------------- Check now if elements are ok ---------------+++++++
+
+ (*_vertices.begin())->onWhat()->model()->writeMSH(meshName);
+
+ double percentage = max_incr * 100.;
+ while(1){
+ std::vector<MElement*> disto;
+ double minD;
+ getDistordedElements(v, 0.5, disto, minD);
+ if (minD < thres){
+ percentage -= 10.;
+ for (std::set<MVertex*>::iterator it = _vertices.begin(); it != _vertices.end(); ++it){
+ double ax, ay, az;
+ myAssembler.getDofValue(*it, 0, _tag, ax);
+ myAssembler.getDofValue(*it, 1, _tag, ay);
+ myAssembler.getDofValue(*it, 2, _tag, az);
+ (*it)->x() -= .1*ax;
+ (*it)->y() -= .1*ay;
+ (*it)->z() -= .1*az;
+ }
+ }
+ else break;
+ }
+
+ delete lsys;
+ return percentage;
+#endif
+ return 0.0;
+}
+
+// uncurve elements that are invalid
+void highOrderTools::ensureMinimumDistorsion (std::vector<MElement*> &all,
+ double threshold){
+ while(1){
+ double minD;
+ std::vector<MElement*> disto;
+ getDistordedElements(all, threshold, disto, minD);
+ if (!disto.size())break;
+ Msg::Info("fixing %d bad curved elements (worst disto %g)",disto.size(),minD);
+ for (int i=0;i<disto.size();i++){
+ ensureMinimumDistorsion(disto[i],threshold);
+ }
+ }
+}
+
+double highOrderTools::applySmoothingTo (std::vector<MElement*> &all,
+ double threshold, bool mixed)
+{
+ int ITER = 0;
+ double minD, FACT = 1.0;
+ std::vector<MElement*> disto;
+ // move the points to their straight sided locations
+ for (unsigned int i = 0; i < all.size(); i++)
+ moveToStraightSidedLocation(all[i]);
+ // apply the displacement
+
+ _gm->writeMSH("straightSided.msh");
+
+
+ double percentage_of_what_is_left = apply_incremental_displacement (1.,all, mixed, -10000000000 ,"sm.msh",all);
+ ensureMinimumDistorsion (all,threshold);
+ return 1.;
+
+ double percentage = 0.0;
+ while(1){
+ char NN[256];
+ sprintf(NN,"smoothing-%d.msh",ITER++);
+ double percentage_of_what_is_left = apply_incremental_displacement (1.,all, mixed, threshold,NN,all);
+ percentage += (1.-percentage) * percentage_of_what_is_left/100.;
+ Msg::Info("The smoother was able to do %3d percent of the motion",(int)(percentage*100.));
+ if (percentage_of_what_is_left == 0.0) break;
+ else if (percentage_of_what_is_left == 100.)break;
+ }
+
+ getDistordedElements(all, 0.3, disto, minD);
+ Msg::Info("iter %d : %d elements / %d distorted min Disto = %g ",ITER,
+ all.size(), disto.size(), minD);
+ return percentage;
+}
+
+extern void printJacobians(GModel *m, const char *nm);
+void highOrderTools::makePosViewWithJacobians (const char *fn){
+ printJacobians(_gm,fn);
+}
+
+#endif
+
diff --git a/Mesh/highOrderTools.h b/Mesh/highOrderTools.h
new file mode 100644
index 0000000000..a740079138
--- /dev/null
+++ b/Mesh/highOrderTools.h
@@ -0,0 +1,92 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _HIGH_ORDER_TOOLS_H_
+#define _HIGH_ORDER_TOOLS_H_
+
+#include <map>
+#include <vector>
+#include "GmshConfig.h"
+#include "GmshMessage.h"
+
+#if defined(HAVE_SOLVER)
+
+#include "SVector3.h"
+#include "fullMatrix.h"
+#include "dofManager.h"
+#include "elasticityTerm.h"
+
+class MVertex;
+class MElement;
+class GFace;
+class GRegion;
+
+class highOrderTools
+{
+ GModel *_gm;
+ const int _tag;
+ // contains the position of vertices in a straight sided mesh
+ std::map<MVertex*,SVector3> _straightSidedLocation;
+ // contains the position of vertices in the best curvilinear mesh
+ // available
+ std::map<MVertex*,SVector3> _targetLocation;
+ int _dim;
+ void _clean()
+ {
+ _straightSidedLocation.clear();
+ _targetLocation.clear();
+ }
+ double smooth_metric_(std::vector<MElement*> & v,
+ GFace *gf,
+ dofManager<double> &myAssembler,
+ std::set<MVertex*> &verticesToMove,
+ elasticityTerm &El);
+ void computeMetricInfo(GFace *gf,
+ MElement *e,
+ fullMatrix<double> &J,
+ fullMatrix<double> &JT,
+ fullVector<double> &D);
+ double apply_incremental_displacement (double max_incr,
+ std::vector<MElement*> & v,
+ bool mixed,
+ double thres,
+ char *meshName,
+ std::vector<MElement*> & disto);
+ void moveToStraightSidedLocation(MElement *) const;
+ void computeStraightSidedPositions ();
+ public:
+ highOrderTools(GModel *gm);
+ highOrderTools(GModel *gm, GModel *mesh, int order);
+ // void applyGlobalSmoothing ();
+ void ensureMinimumDistorsion(MElement *e, double threshold);
+ void ensureMinimumDistorsion (std::vector<MElement*> &all, double threshold);
+ void ensureMinimumDistorsion (double threshold);
+ double applySmoothingTo (GFace *gf, double tres = 0.1, bool mixed = false);
+ void applySmoothingTo(std::vector<MElement*> & all, GFace *gf);
+ double applySmoothingTo (std::vector<MElement*> &all, double threshold, bool mixed);
+ inline SVector3 getSSL(MVertex *v) const
+ {
+ std::map<MVertex*,SVector3>::const_iterator it = _straightSidedLocation.find(v);
+ if (it != _straightSidedLocation.end())return it->second;
+ else return SVector3(v->x(),v->y(),v->z());
+ }
+ void makePosViewWithJacobians(const char *nm);
+};
+
+#else
+
+class highOrderTools
+{
+ public:
+ highOrderTools(GModel *gm);
+ {
+ Msg::Error("Gmsh has to be compiled with solver support to use highOrderSmoother");
+ }
+ void applyGlobalSmoothing (){}
+};
+
+#endif
+
+#endif
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index bf064ae572..3705a67305 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -231,6 +231,12 @@ void buildEdgeToTriangle(std::vector<MTriangle*> &tris, e2t_cont &adj)
buildEdgeToElement(tris, adj);
}
+void buildEdgeToElements(std::vector<MElement*> &tris, e2t_cont &adj)
+{
+ adj.clear();
+ buildEdgeToElement(tris, adj);
+}
+
void buildListOfEdgeAngle(e2t_cont adj, std::vector<edge_angle> &edges_detected,
std::vector<edge_angle> &edges_lonly)
{
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index 98d6b6f17e..988d82a96a 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -54,6 +54,8 @@ void buildVertexToTriangle(std::vector<MTriangle*> &, v2t_cont &adj);
void buildEdgeToTriangle(std::vector<MTriangle*> &, e2t_cont &adj);
void buildListOfEdgeAngle(e2t_cont adj, std::vector<edge_angle> &edges_detected,
std::vector<edge_angle> &edges_lonly);
+void buildEdgeToElements(std::vector<MElement*> &tris, e2t_cont &adj);
+
void laplaceSmoothing(GFace *gf, int niter=1);
void edgeSwappingLawson(GFace *gf);
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index bfa67191ca..e4bed315d4 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -212,20 +212,22 @@ double mesh_functional_distorsion(MElement *t, double u, double v)
double mesh_functional_distorsion_p2(MTriangle *t)
{
- double d = mesh_functional_distorsion(t,0.,0.);
- d = std::min(d,mesh_functional_distorsion(t,1.,0.));
- d = std::min(d,mesh_functional_distorsion(t,0.,1.));
- d = std::min(d,mesh_functional_distorsion(t,.5,0.));
- d = std::min(d,mesh_functional_distorsion(t,0.,0.5));
- d = std::min(d,mesh_functional_distorsion(t,.5,0.5));
- return d;
+ double d1 =mesh_functional_distorsion(t,0.0,0.0);
+ double d2 =mesh_functional_distorsion(t,1.0,0.0);
+ double d3 =mesh_functional_distorsion(t,0.0,1.0);
+ double d4 =mesh_functional_distorsion(t,0.5,0.0);
+ double d5 =mesh_functional_distorsion(t,0.5,0.5);
+ double d6 =mesh_functional_distorsion(t,0.0,0.5);
+ double d[6] = {d1,d2,d3,4*d4-d1-d2,4*d5-d2-d3,4*d6-d1-d3};
+
+ return *std::min_element(d,d+6);
}
double qmDistorsionOfMapping (MTriangle *e)
{
// return 1.0;
if (e->getPolynomialOrder() == 1) return 1.0;
- // if (e->getPolynomialOrder() == 2) return mesh_functional_distorsion_p2(e);
+ if (e->getPolynomialOrder() == 2) return mesh_functional_distorsion_p2(e);
IntPt *pts;
int npts;
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index e732bcd9b0..8071ad32d1 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -835,7 +835,7 @@ DocRecord::~DocRecord()
for(int i = 0; i < numPoints; i++)
delete [] _adjacencies[i].t;
delete _adjacencies;
- }
+ }
}
void DocRecord::MakeMeshWithPoints()
@@ -853,6 +853,59 @@ void DocRecord::Voronoi()
ConvertDListToVoronoiData();
}
+void DocRecord::recur_tag_triangles (int iTriangle,
+ std::set<int> &taggedTriangles,
+ std::map< std::pair<void*,void*>, std::pair<int,int> > & edgesToTriangles){
+
+ if (taggedTriangles.find(iTriangle) != taggedTriangles.end())return;
+
+ taggedTriangles.insert(iTriangle);
+
+ int a = triangles[iTriangle].a;
+ int b = triangles[iTriangle].b;
+ int c = triangles[iTriangle].c;
+ PointRecord *p[3] = {&points[a],&points[b],&points[c]};
+
+ for (int j=0;j<3;j++){
+ if (!lookForBoundaryEdge(p[j]->data,p[(j+1)%3]->data)){
+ std::pair<void*,void*> ab = std::make_pair(std::min(p[j]->data,p[(j+1)%3]->data),std::max(p[j]->data,p[(j+1)%3]->data));
+ std::pair<int,int> & adj = edgesToTriangles[ab];
+ if (iTriangle == adj.first && adj.second != -1)recur_tag_triangles(adj.second, taggedTriangles, edgesToTriangles);
+ else recur_tag_triangles(adj.first, taggedTriangles, edgesToTriangles);
+ }
+ }
+}
+
+
+std::set<int> DocRecord::tagInterior (double x, double y){
+ std::map< std::pair<void*,void*>, std::pair<int,int> > edgesToTriangles;
+ int iFirst = 1;
+ for (int i=0;i<numTriangles;i++){
+ int a = triangles[i].a;
+ int b = triangles[i].b;
+ int c = triangles[i].c;
+ PointRecord *p[3] = {&points[a],&points[b],&points[c]};
+
+ // TODO : regarde quel triangle contient x y et retiens le dans iFirst;
+
+ for (int j=0;j<3;j++){
+ std::pair<void*,void*> ab = std::make_pair(std::min(p[j]->data,p[(j+1)%3]->data),std::max(p[j]->data,p[(j+1)%3]->data));
+ std::map< std::pair<void*,void*>, std::pair<int,int> > :: iterator it = edgesToTriangles.find(ab);
+ if (it == edgesToTriangles.end()){
+ edgesToTriangles[ab] = std::make_pair(i,-1);
+ }
+ else {
+ it->second.second = i;
+ }
+ }
+ }
+
+ std::set<int> taggedTriangles;
+ recur_tag_triangles (iFirst, taggedTriangles, edgesToTriangles);
+ return taggedTriangles;
+}
+
+
void DocRecord::setPoints(fullMatrix<double> *p)
{
if (numPoints != p->size1())throw;
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index 98c6b24622..a04b72cbb2 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -89,11 +89,26 @@ class DocRecord{
PointRecord *points; // points to triangulate
int numTriangles; // number of triangles
Triangle *triangles; // 2D results
+ std::set<std::pair<void*,void*> > boundaryEdges;
DocRecord(int n);
double &x(int i){ return points[i].where.h; }
double &y(int i){ return points[i].where.v; }
void* &data(int i){ return points[i].data; }
- void setPoints(fullMatrix<double> *p);
+ void addBoundaryEdge (void* p1, void* p2){
+ void *a = (p1 < p2) ? p1 : p2;
+ void *b = (p1 > p2) ? p1 : p2;
+ boundaryEdges.insert(std::make_pair(a,b));
+ }
+ bool lookForBoundaryEdge(void *p1, void* p2){
+ void *a = (p1 < p2) ? p1 : p2;
+ void *b = (p1 > p2) ? p1 : p2;
+ std::set<std::pair<void*,void*> >::iterator it = boundaryEdges.find(std::make_pair(a,b));
+ return it != boundaryEdges.end();
+ }
+ void recur_tag_triangles (int iTriangle,
+ std::set<int> &taggedTriangles,
+ std::map< std::pair<void*,void*>, std::pair<int,int> > & edgesToTriangles);
+ void setPoints(fullMatrix<double> *p);
~DocRecord();
void MakeMeshWithPoints();
void Voronoi ();
@@ -108,6 +123,7 @@ class DocRecord{
void remove_all();
void add_point(double,double,GFace*);
PointNumero *ConvertDlistToArray(DListPeek *dlist, int *n);
+ std::set<int> tagInterior (double x, double y);
};
void centroidOfOrientedBox(std::vector<SPoint2> &pts,
diff --git a/Solver/elasticityTerm.cpp b/Solver/elasticityTerm.cpp
index 945c7c5f80..dc034b2a29 100644
--- a/Solver/elasticityTerm.cpp
+++ b/Solver/elasticityTerm.cpp
@@ -143,7 +143,7 @@ void elasticityMixedTerm::elementMatrix(SElement *se, fullMatrix<double> &m) con
{
setPolynomialBasis(se);
MElement *e = se->getMeshElement();
- int integrationOrder = 3 * (e->getPolynomialOrder()) + 2;
+ int integrationOrder = 4 * (e->getPolynomialOrder()) + 2;
int npts;
IntPt *GP;
double jac[3][3];
@@ -158,6 +158,7 @@ void elasticityMixedTerm::elementMatrix(SElement *se, fullMatrix<double> &m) con
fullMatrix<double> KUG(3 * _sizeN,_sizeM);
fullMatrix<double> KPG(_sizeM,_sizeM);
fullMatrix<double> KGG(_sizeM,_sizeM);
+ fullMatrix<double> KPP(_sizeM,_sizeM); // stabilization
double FACT = _E / ((1 + _nu)*(1.-2*_nu));
double C11 = FACT * (1 - _nu) ;
@@ -192,7 +193,7 @@ void elasticityMixedTerm::elementMatrix(SElement *se, fullMatrix<double> &m) con
for (int j = 0; j < 6; j++)
H(i, j) = C[i][j];
- double _dimension = 3.0;
+ double _dimension = 2.0;
for (int i = 0; i < npts; i++){
const double u = GP[i].pt[0];
@@ -206,33 +207,69 @@ void elasticityMixedTerm::elementMatrix(SElement *se, fullMatrix<double> &m) con
B.setAll(0.);
BT.setAll(0.);
+
+ const double K = .0000002*weight * detJ * pow (detJ,2/_dimension); // the second detJ is for stabilization
+
+ for (int j = 0; j < _sizeM; j++){
+ for (int k = 0; k < _sizeM; k++){
+ KPP(j, k) += (Grads[j][0] * Grads[k][0] +
+ Grads[j][1] * Grads[k][1] +
+ Grads[j][2] * Grads[k][2])* K;
+ }
+ }
+
+ const double K2 = weight * detJ;
+
+ for (int j = 0; j < _sizeN; j++){
+ for (int k = 0; k < _sizeM; k++){
+ KUP(j+0*_sizeN, k) += Grads[j][0] * SF[k] * K2;
+ KUP(j+1*_sizeN, k) += Grads[j][1] * SF[k] * K2;
+ KUP(j+2*_sizeN, k) += Grads[j][2] * SF[k] * K2;
+ }
+ }
+
+ const double K3 = weight * detJ * _E/(2*(1+_nu));
+ /*
+ for (int j = 0; j < _sizeN; j++){
+ for (int k = 0; k < _sizeN; k++){
+ const double fa = (Grads[j][0] * Grads[k][0] +
+ Grads[j][1] * Grads[k][1] +
+ Grads[j][2] * Grads[k][2]) * K3;
+ KUU(j+0*_sizeN, k+0*_sizeN) += fa;
+ KUU(j+1*_sizeN, k+1*_sizeN) += fa;
+ KUU(j+2*_sizeN, k+2*_sizeN) += fa;
+ }
+ }
+ */
+
+
for (int j = 0; j < _sizeN; j++){
// printf("%g %g %g\n",Grads[j][0],Grads[j][1],Grads[j][2]);
BDILT(j, 0) = BDIL(0, j) = Grads[j][0]/_dimension;
- BT(j, 0) = B(0, j) = Grads[j][0]-Grads[j][0]/_dimension;
- BT(j, 1) = B(1, j) = -Grads[j][1]/_dimension;
- BT(j, 2) = B(2, j) = -Grads[j][2]/_dimension;
+ BT(j, 0) = B(0, j) = Grads[j][0];//-Grads[j][0]/_dimension;
+ // BT(j, 1) = B(1, j) = -Grads[j][1]/_dimension;
+ // BT(j, 2) = B(2, j) = -Grads[j][2]/_dimension;
BT(j, 3) = B(3, j) = Grads[j][1];
BT(j, 4) = B(4, j) = Grads[j][2];
BDILT(j + _sizeN, 0) = BDIL(0, j + _sizeN) = Grads[j][1]/_dimension;
- BT(j + _sizeN, 0) = B(0, j + _sizeN) = -Grads[j][0]/_dimension;
- BT(j + _sizeN, 1) = B(1, j + _sizeN) = Grads[j][1]-Grads[j][1]/_dimension;
- BT(j + _sizeN, 2) = B(2, j + _sizeN) = -Grads[j][2]/_dimension;
+ // BT(j + _sizeN, 0) = B(0, j + _sizeN) = -Grads[j][0]/_dimension;
+ BT(j + _sizeN, 1) = B(1, j + _sizeN) = Grads[j][1];//-Grads[j][1]/_dimension;
+ // BT(j + _sizeN, 2) = B(2, j + _sizeN) = -Grads[j][2]/_dimension;
BT(j + _sizeN, 3) = B(3, j + _sizeN) = Grads[j][0];
BT(j + _sizeN, 5) = B(5, j + _sizeN) = Grads[j][2];
BDILT(j + 2 * _sizeN, 0) = BDIL(0, j + 2 * _sizeN) = Grads[j][2]/_dimension;
- BT(j + 2 * _sizeN, 0) = B(0, j + 2 * _sizeN) = -Grads[j][0]/_dimension;
- BT(j + 2 * _sizeN, 1) = B(1, j + 2 * _sizeN) = -Grads[j][1]/_dimension;
- BT(j + 2 * _sizeN, 2) = B(2, j + 2 * _sizeN) = Grads[j][2]-Grads[j][2]/_dimension;
+ // BT(j + 2 * _sizeN, 0) = B(0, j + 2 * _sizeN) = -Grads[j][0]/_dimension;
+ // BT(j + 2 * _sizeN, 1) = B(1, j + 2 * _sizeN) = -Grads[j][1]/_dimension;
+ BT(j + 2 * _sizeN, 2) = B(2, j + 2 * _sizeN) = Grads[j][2];//-Grads[j][2]/_dimension;
BT(j + 2 * _sizeN, 4) = B(4, j + 2 * _sizeN) = Grads[j][1];
BT(j + 2 * _sizeN, 5) = B(5, j + 2 * _sizeN) = Grads[j][0];
@@ -252,17 +289,17 @@ void elasticityMixedTerm::elementMatrix(SElement *se, fullMatrix<double> &m) con
}
KUU.gemm(BTH, B, weight * detJ);
- KUP.gemm(BDILT, M, weight * detJ);
- KPG.gemm(MT, M, -weight * detJ);
- KUG.gemm(trBTH, M, weight * detJ/_dimension);
- KGG.gemm(MT, M, weight * detJ * (_K)/(_dimension*_dimension));
+ // KUP.gemm(BDILT, M, weight * detJ);
+ // KPG.gemm(MT, M, -weight * detJ);
+ // KUG.gemm(trBTH, M, weight * detJ/_dimension);
+ // KGG.gemm(MT, M, weight * detJ * (_K)/(_dimension*_dimension));
}
/*
3*_sizeN _sizeM _sizeM
KUU KUP KUG
- m = KUP^t 0 KPG
+ m = KUP^t KPP KPG
KUG^T KPG^t KGG
*/
@@ -288,6 +325,8 @@ void elasticityMixedTerm::elementMatrix(SElement *se, fullMatrix<double> &m) con
m(i+3*_sizeN+_sizeM,j+3*_sizeN+_sizeM) = KGG(i,j); // assemble KGG
+ m(i+3*_sizeN,j+3*_sizeN) = KPP(i,j); // assemble KPP
+
}
}
// m.print("Mixed");
--
GitLab