From 2bd27f429ca134bd4f092083f8fbcf5ad53d8a99 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 23 Jan 2011 15:48:26 +0000
Subject: [PATCH] basic 2D boundary layers (using old code)
---
Mesh/BoundaryLayers.cpp | 104 +++++++++++++++++++++++++++++-----------
Mesh/meshGEdge.cpp | 89 +++++-----------------------------
2 files changed, 86 insertions(+), 107 deletions(-)
diff --git a/Mesh/BoundaryLayers.cpp b/Mesh/BoundaryLayers.cpp
index 760a7b3c91..9320f81680 100644
--- a/Mesh/BoundaryLayers.cpp
+++ b/Mesh/BoundaryLayers.cpp
@@ -4,6 +4,7 @@
// bugs and problems to <gmsh@geuz.org>.
#include "GModel.h"
+#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "BoundaryLayers.h"
@@ -38,16 +39,17 @@ static void addExtrudeNormals(std::vector<T*> &elements, int invert,
else{ // get extrusion data from Gouraud-shaded element normals
for(unsigned int i = 0; i < elements.size(); i++){
MElement *ele = elements[i];
- for(int j = 0; j < ele->getNumFaces(); j++){
- MFace fac = ele->getFace(j);
- SVector3 n = fac.normal();
- if(invert) n *= -1.;
- if(n[0] || n[1] || n[2]){
- double nn[3] = {n[0], n[1], n[2]};
- for(int k = 0; k < fac.getNumVertices(); k++){
- MVertex *v = fac.getVertex(k);
- ExtrudeParams::normals->add(v->x(), v->y(), v->z(), 3, nn);
- }
+ SVector3 n(0, 0, 0);
+ if(ele->getDim() == 2)
+ n = ele->getFace(0).normal();
+ else if(ele->getDim() == 1) // FIXME only valid in XY-plane
+ n = crossprod(ele->getEdge(0).tangent(), SVector3(0, 0, 1));
+ if(invert) n *= -1.;
+ if(n[0] || n[1] || n[2]){
+ double nn[3] = {n[0], n[1], n[2]};
+ for(int k = 0; k < ele->getNumVertices(); k++){
+ MVertex *v = ele->getVertex(k);
+ ExtrudeParams::normals->add(v->x(), v->y(), v->z(), 3, nn);
}
}
}
@@ -58,8 +60,29 @@ int Mesh2DWithBoundaryLayers(GModel *m)
{
std::set<GFace*> sourceFaces, otherFaces;
std::set<GEdge*> sourceEdges, otherEdges;
- std::map<int, bool> sourceFaceInvert;
- std::map<int, int> sourceUseView;
+ std::map<int, bool> sourceFaceInvert, sourceEdgeInvert;
+ std::map<int, int> sourceFaceUseView, sourceEdgeUseView;
+
+ // 2D boundary layers
+ for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++){
+ GEdge *ge = *it;
+ if(ge->getNativeType() == GEntity::GmshModel &&
+ ge->geomType() == GEntity::BoundaryLayerCurve){
+ ExtrudeParams *ep = ge->meshAttributes.extrude;
+ if(ep && ep->mesh.ExtrudeMesh && ep->geo.Mode == COPIED_ENTITY){
+ GEdge *from = m->getEdgeByTag(std::abs(ep->geo.Source));
+ if(!from){
+ Msg::Error("Unknown source curve %d for boundary layer", ep->geo.Source);
+ return 0;
+ }
+ if(ep->geo.Source < 0) sourceEdgeInvert[from->tag()] = true;
+ if(ep->mesh.ViewIndex >= 0) sourceEdgeUseView[from->tag()] = ep->mesh.ViewIndex;
+ sourceEdges.insert(from);
+ }
+ }
+ }
+
+ // 3D boundary layers
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); it++){
GFace *gf = *it;
if(gf->getNativeType() == GEntity::GmshModel &&
@@ -72,14 +95,16 @@ int Mesh2DWithBoundaryLayers(GModel *m)
return 0;
}
if(ep->geo.Source < 0) sourceFaceInvert[from->tag()] = true;
- if(ep->mesh.ViewIndex >= 0) sourceUseView[from->tag()] = ep->mesh.ViewIndex;
+ if(ep->mesh.ViewIndex >= 0) sourceFaceUseView[from->tag()] = ep->mesh.ViewIndex;
sourceFaces.insert(from);
std::list<GEdge*> e = from->edges();
sourceEdges.insert(e.begin(), e.end());
}
}
}
- if(sourceFaces.empty()) return 0;
+
+ if(sourceEdges.empty() && sourceFaces.empty())
+ return 0;
// compute set of non-source edges and faces
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); it++)
@@ -96,27 +121,48 @@ int Mesh2DWithBoundaryLayers(GModel *m)
// extrusion)
std::for_each(sourceFaces.begin(), sourceFaces.end(), orientMeshGFace());
- // compute a normal field for all the source faces
+ // compute a normal field for all the source edges or faces
if(ExtrudeParams::normals) delete ExtrudeParams::normals;
ExtrudeParams::normals = new smooth_data();
- for(std::set<GFace*>::iterator it = sourceFaces.begin();
- it != sourceFaces.end(); it++){
- GFace *gf = *it;
- int invert = sourceFaceInvert.count(gf->tag());
- OctreePost *octree = 0;
+
+ if(sourceFaces.empty()){
+ for(std::set<GEdge*>::iterator it = sourceEdges.begin();
+ it != sourceEdges.end(); it++){
+ GEdge *ge = *it;
+ int invert = sourceEdgeInvert.count(ge->tag());
+ OctreePost *octree = 0;
#if defined(HAVE_POST)
- if(sourceUseView.count(gf->tag())){
- int index = sourceUseView[gf->tag()];
- if(index >= 0 && index < PView::list.size())
- octree = new OctreePost(PView::list[index]);
- else
- Msg::Error("Unknown View[%d]: using normals instead", index);
+ if(sourceEdgeUseView.count(ge->tag())){
+ int index = sourceEdgeUseView[ge->tag()];
+ if(index >= 0 && index < PView::list.size())
+ octree = new OctreePost(PView::list[index]);
+ else
+ Msg::Error("Unknown View[%d]: using normals instead", index);
+ }
+#endif
+ addExtrudeNormals(ge->lines, invert, octree);
}
+ }
+ else{
+ for(std::set<GFace*>::iterator it = sourceFaces.begin();
+ it != sourceFaces.end(); it++){
+ GFace *gf = *it;
+ int invert = sourceFaceInvert.count(gf->tag());
+ OctreePost *octree = 0;
+#if defined(HAVE_POST)
+ if(sourceFaceUseView.count(gf->tag())){
+ int index = sourceFaceUseView[gf->tag()];
+ if(index >= 0 && index < PView::list.size())
+ octree = new OctreePost(PView::list[index]);
+ else
+ Msg::Error("Unknown View[%d]: using normals instead", index);
+ }
#endif
- addExtrudeNormals(gf->triangles, invert, octree);
- addExtrudeNormals(gf->quadrangles, invert, octree);
+ addExtrudeNormals(gf->triangles, invert, octree);
+ addExtrudeNormals(gf->quadrangles, invert, octree);
+ }
}
- if(sourceUseView.empty())
+ if(sourceEdgeUseView.empty() && sourceFaceUseView.empty())
ExtrudeParams::normals->normalize();
// set the position of boundary layer points using the smooth normal
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 89ea5c5e04..0d6f918420 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -20,66 +20,6 @@ typedef struct {
double t, lc, p;
} IntPoint;
-struct xi2lc {
- double xi, lc;
- xi2lc(const double &_xi, const double _lc)
- : xi(_xi), lc(_lc)
- {
- }
- bool operator < (const xi2lc &other) const
- {
- return xi < other.xi;
- }
-};
-
-static std::vector<xi2lc> interpLc;
-
-static void buildInterpLc(const std::vector<IntPoint> &lcPoints)
-{
- IntPoint p;
- interpLc.clear();
- for(unsigned int i = 0; i < lcPoints.size(); i++){
- p = lcPoints[i];
- interpLc.push_back(xi2lc(p.t, p.lc));
- }
-}
-
-static double F_Lc_usingInterpLc(GEdge *ge, double t)
-{
- std::vector<xi2lc>::iterator it = std::lower_bound(interpLc.begin(),
- interpLc.end(), xi2lc(t, 0));
- double t1 = it->xi;
- double l1 = it->lc;
- it++;
- if(it == interpLc.end()) return l1;
- double t2 = it->xi;
- double l2 = it->lc;
- double l = l1 + ((t - t1) / (t2 - t1)) * (l2 - l1);
- return l;
-}
-
-static double F_Lc_usingInterpLcBis(GEdge *ge, double t)
-{
- GPoint p = ge->point(t);
- double lc_here;
-
- Range<double> bounds = ge->parBounds(0);
- double t_begin = bounds.low();
- double t_end = bounds.high();
-
- SVector3 der = ge->firstDer(t);
- const double d = norm(der);
-
- if(t == t_begin)
- lc_here = BGM_MeshSize(ge->getBeginVertex(), t, 0, p.x(), p.y(), p.z());
- else if(t == t_end)
- lc_here = BGM_MeshSize(ge->getEndVertex(), t, 0, p.x(), p.y(), p.z());
- else
- lc_here = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
-
- return d / lc_here;
-}
-
static double F_Lc(GEdge *ge, double t)
{
GPoint p = ge->point(t);
@@ -122,7 +62,7 @@ static double F_Lc_aniso(GEdge *ge, double t)
SVector3 der = ge->firstDer(t);
- double lSquared = dot (der,lc_here,der);
+ double lSquared = dot(der, lc_here, der);
// der.normalize();
// printf("in the function %g n %g %g\n", sqrt(lSquared),der.x(),der.y());
@@ -367,25 +307,18 @@ void meshGEdge::operator() (GEdge *ge)
N = ge->meshAttributes.nbPointsTransfinite;
}
else{
- /*if(CTX::instance()->mesh.lcIntegrationPrecision > 1.e-2){ JF says this code was a test but it is not useful
- std::vector<IntPoint> lcPoints;
- Integration(ge, t_begin, t_end, F_Lc_usingInterpLcBis, lcPoints,
- CTX::instance()->mesh.lcIntegrationPrecision);
- buildInterpLc(lcPoints);
- a = Integration(ge, t_begin, t_end, F_Lc_usingInterpLc, Points, 1.e-8);
- }
- else*/{
- if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG)
- a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
- CTX::instance()->mesh.lcIntegrationPrecision);
- else
- a = Integration(ge, t_begin, t_end, F_Lc, Points,
- CTX::instance()->mesh.lcIntegrationPrecision);
- }
+ if (CTX::instance()->mesh.algo2d == ALGO_2D_BAMG)
+ a = Integration(ge, t_begin, t_end, F_Lc_aniso, Points,
+ CTX::instance()->mesh.lcIntegrationPrecision);
+ else
+ a = Integration(ge, t_begin, t_end, F_Lc, Points,
+ CTX::instance()->mesh.lcIntegrationPrecision);
N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
}
-
- if(CTX::instance()->mesh.algoRecombine == 1 && N%2 == 0) N++;
+
+ // force odd number of points for if blossom is used for
+ // recombination
+ if(CTX::instance()->mesh.algoRecombine == 1 && N % 2 == 0) N++;
// printFandPrimitive(ge->tag(),Points);
--
GitLab