Select Git revision
meshGEdge.cpp 23.96 KiB
// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
#include "GmshConfig.h"
#include "GModel.h"
#include "discreteEdge.h"
#include "meshGEdge.h"
#include "GEdge.h"
#include "MLine.h"
#include "BackgroundMeshTools.h"
#include "boundaryLayersData.h"
#include "Numeric.h"
#include "GmshMessage.h"
#include "Context.h"
#include "STensor3.h"
#include "Field.h"
#include "OS.h"
typedef struct {
int Num;
// t is the local coordinate of the point
// lc is x'(t)/h(x(t))
// p is the value of the primitive
// xp is the norm of the tangent vector
double t, lc, p, xp;
} IntPoint;
static double smoothPrimitive(GEdge *ge, double alpha,
std::vector<IntPoint> &Points)
{
int ITER = 0;
while(1) {
int count1 = 0;
int count2 = 0;
// use a gauss-seidel iteration; iterate forward and then backward;
// convergence is usually very fast
for(std::size_t i = 1; i < Points.size(); i++) {
double dh =
(Points[i].xp / Points[i].lc - Points[i - 1].xp / Points[i - 1].lc);
double dt = Points[i].t - Points[i - 1].t;
double dhdt = dh / dt;
if(dhdt / Points[i].xp > (alpha - 1.) * 1.01) {
double hnew =
Points[i - 1].xp / Points[i - 1].lc + dt * (alpha - 1) * Points[i].xp;
Points[i].lc = Points[i].xp / hnew;
count1++;
}
}
for(int i = Points.size() - 1; i > 0; i--) {
double dh =
(Points[i - 1].xp / Points[i - 1].lc - Points[i].xp / Points[i].lc);
double dt = std::abs(Points[i - 1].t - Points[i].t);
double dhdt = dh / dt;
if(dhdt / Points[i - 1].xp > (alpha - 1.) * 1.01) {
double hnew =
Points[i].xp / Points[i].lc + dt * (alpha - 1) * Points[i].xp;
Points[i - 1].lc = Points[i - 1].xp / hnew;
count2++;
}
}
++ITER;
if(ITER > 2000) break;
if(!(count1 + count2)) break;
}
// recompute the primitive
for(std::size_t i = 1; i < Points.size(); i++) {
IntPoint &pt2 = Points[i];
IntPoint &pt1 = Points[i - 1];
pt2.p = pt1.p + (pt2.t - pt1.t) * 0.5 * (pt2.lc + pt1.lc);
}
return Points[Points.size() - 1].p;
}
struct F_LcB {
double operator()(GEdge *ge, double t)
{
GPoint p = ge->point(t);
return BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
}
};
struct F_Lc {
double operator()(GEdge *ge, double t)
{
GPoint p = ge->point(t);
Range<double> bounds = ge->parBounds(0);
double t_begin = bounds.low();
double t_end = bounds.high();
double lc_here;
if(t == t_begin && ge->getBeginVertex())
lc_here = BGM_MeshSize(ge->getBeginVertex(), t, 0, p.x(), p.y(), p.z());
else if(t == t_end && ge->getEndVertex())
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());
SVector3 der = ge->firstDer(t);
return norm(der) / lc_here;
}
};
struct F_Lc_aniso {
double operator()(GEdge *ge, double t)
{
GPoint p = ge->point(t);
SMetric3 lc_here;
Range<double> bounds = ge->parBounds(0);
double t_begin = bounds.low();
double t_end = bounds.high();
if(t == t_begin && ge->getBeginVertex())
lc_here = BGM_MeshMetric(ge->getBeginVertex(), t, 0, p.x(), p.y(), p.z());
else if(t == t_end && ge->getEndVertex())
lc_here = BGM_MeshMetric(ge->getEndVertex(), t, 0, p.x(), p.y(), p.z());
else
lc_here = BGM_MeshMetric(ge, t, 0, p.x(), p.y(), p.z());
FieldManager *fields = ge->model()->getFields();
for(int i = 0; i < fields->getNumBoundaryLayerFields(); ++i) {
Field *bl_field = fields->get(fields->getBoundaryLayerField(i));
if(bl_field == NULL) continue;
BoundaryLayerField *blf = dynamic_cast<BoundaryLayerField *>(bl_field);
if(blf->isEdgeBL(ge->tag())) break;
SMetric3 lc_bgm;
blf->computeFor1dMesh(p.x(), p.y(), p.z(), lc_bgm);
lc_here = intersection_conserveM1(lc_here, lc_bgm);
}
SVector3 der = ge->firstDer(t);
return std::sqrt(dot(der, lc_here, der));
}
};
struct F_Transfinite {
double operator()(GEdge *ge, double t_)
{
double length = ge->length();
if(length == 0.0) {
Msg::Error("Zero-length curve %d in transfinite mesh", ge->tag());
return 1.;
}
SVector3 der = ge->firstDer(t_);
double d = norm(der);
double coef = ge->meshAttributes.coeffTransfinite;
int type = ge->meshAttributes.typeTransfinite;
int nbpt = ge->meshAttributes.nbPointsTransfinite;
if(CTX::instance()->mesh.flexibleTransfinite &&
CTX::instance()->mesh.lcFactor)
nbpt /= CTX::instance()->mesh.lcFactor;
Range<double> bounds = ge->parBounds(0);
double t_begin = bounds.low();
double t_end = bounds.high();
double t = (t_ - t_begin) / (t_end - t_begin);
double val;
if(coef <= 0.0 || coef == 1.0) {
// coef < 0 should never happen
val = d * coef / ge->length();
}
else {
switch(std::abs(type)) {
case 1: // Geometric progression ar^i; Sum of n terms = length = a
// (r^n-1)/(r-1)
{
double r = (gmsh_sign(type) >= 0) ? coef : 1. / coef;
double a = length * (r - 1.) / (std::pow(r, nbpt - 1.) - 1.);
int i = (int)(std::log(t * length / a * (r - 1.) + 1.) / std::log(r));
val = d / (a * std::pow(r, (double)i));
} break;
case 2: // Bump
{
double a;
if(coef > 1.0) {
a = -4. * std::sqrt(coef - 1.) *
std::atan2(1.0, std::sqrt(coef - 1.)) / ((double)nbpt * length);
}
else {
a = 2. * std::sqrt(1. - coef) *
std::log(std::abs((1. + 1. / std::sqrt(1. - coef)) /
(1. - 1. / std::sqrt(1. - coef)))) /
((double)nbpt * length);
}
double b = -a * length * length / (4. * (coef - 1.));
val = d / (-a * std::pow(t * length - (length)*0.5, 2) + b);
break;
}
default:
Msg::Warning("Unknown case in Transfinite Line mesh");
val = 1.;
break;
}
}
return val;
}
};
struct F_One {
double operator()(GEdge *ge, double t)
{
SVector3 der = ge->firstDer(t);
return norm(der);
}
};
static double trapezoidal(IntPoint *const P1, IntPoint *const P2)
{
return 0.5 * (P1->lc + P2->lc) * (P2->t - P1->t);
}
template <typename function>
static void RecursiveIntegration(GEdge *ge, IntPoint *from, IntPoint *to,
function f, std::vector<IntPoint> &Points,
double Prec, int *depth)
{
IntPoint P, p1;
(*depth)++;
P.t = 0.5 * (from->t + to->t);
P.lc = f(ge, P.t);
double const val1 = trapezoidal(from, to);
double const val2 = trapezoidal(from, &P);
double const val3 = trapezoidal(&P, to);
double const err = std::abs(val1 - val2 - val3);
if(((err < Prec) && (*depth > 6)) || (*depth > 25)) {
p1 = Points.back();
P.p = p1.p + val2;
Points.push_back(P);
p1 = Points.back();
to->p = p1.p + val3;
Points.push_back(*to);
}
else {
RecursiveIntegration(ge, from, &P, f, Points, Prec, depth);
RecursiveIntegration(ge, &P, to, f, Points, Prec, depth);
}
(*depth)--;
}
template <typename function>
static double Integration(GEdge *ge, double t1, double t2, function f,
std::vector<IntPoint> &Points, double Prec)
{
IntPoint from, to;
int depth = 0;
from.t = t1;
from.lc = f(ge, from.t);
from.p = 0.0;
Points.push_back(from);
to.t = t2;
to.lc = f(ge, to.t);
RecursiveIntegration(ge, &from, &to, f, Points, Prec, &depth);
return Points.back().p;
}
void copyMesh(GEdge *from, GEdge *to, int direction)
{
if(!from->getBeginVertex() || !from->getEndVertex() ||
!to->getBeginVertex() || !to->getEndVertex()){
Msg::Error("Cannot copy mesh on curves without begin/end points");
return;
}
Range<double> u_bounds = from->parBounds(0);
double u_min = u_bounds.low();
double u_max = u_bounds.high();
Range<double> to_u_bounds = to->parBounds(0);
double to_u_min = to_u_bounds.low();
// include begin and end point to avoid conflicts when realigning
MVertex *vt0 = to->getBeginVertex()->mesh_vertices[0];
MVertex *vt1 = to->getEndVertex()->mesh_vertices[0];
MVertex *vs0 = from->getBeginVertex()->mesh_vertices[0];
MVertex *vs1 = from->getEndVertex()->mesh_vertices[0];
to->correspondingVertices[vt0] = direction > 0 ? vs0 : vs1;
to->correspondingVertices[vt1] = direction > 0 ? vs1 : vs0;
for(std::size_t i = 0; i < from->mesh_vertices.size(); i++) {
int index = (direction < 0) ? (from->mesh_vertices.size() - 1 - i) : i;
MVertex *v = from->mesh_vertices[index];
double u;
v->getParameter(0, u);
double newu =
(direction > 0) ? (u - u_min + to_u_min) : (u_max - u + to_u_min);
GPoint gp = to->point(newu);
MEdgeVertex *vv = new MEdgeVertex(gp.x(), gp.y(), gp.z(), to, newu);
to->mesh_vertices.push_back(vv);
to->correspondingVertices[vv] = v;
}
for(std::size_t i = 0; i < to->mesh_vertices.size() + 1; i++) {
MVertex *v0 = (i == 0) ? to->getBeginVertex()->mesh_vertices[0] :
to->mesh_vertices[i - 1];
MVertex *v1 = (i == to->mesh_vertices.size()) ?
to->getEndVertex()->mesh_vertices[0] :
to->mesh_vertices[i];
to->lines.push_back(new MLine(v0, v1));
}
}
void deMeshGEdge::operator()(GEdge *ge)
{
if(ge->geomType() == GEntity::DiscreteCurve){
if(!static_cast<discreteEdge *>(ge)->haveParametrization())
return;
}
ge->deleteMesh();
ge->meshStatistics.status = GEdge::PENDING;
ge->correspondingVertices.clear();
ge->correspondingHOPoints.clear();
}
/*
static void printFandPrimitive(int tag, std::vector<IntPoint> &Points)
{
char name[256];
sprintf(name, "line%d.dat", tag);
FILE *f = Fopen(name, "w");
if(!f) return;
double l = 0;
for (std::size_t i = 0; i < Points.size(); i++){
const IntPoint &P = Points[i];
if(i) l += (P.t - Points[i-1].t) * P.xp;
fprintf(f, "%g %g %g %g %g\n", P.t, P.xp/P.lc, P.p, P.lc, l);
}
fclose(f);
}
*/
// new algo for recombining + splitting
static int increaseN(int N)
{
if(((N + 1) / 2 - 1) % 2 != 0) return N + 2;
return N;
}
// ensure not to have points that are too close to each other.
// can be caused by a coarse 1D mesh or by a noisy curve
static void filterPoints(GEdge *ge, int nMinimumPoints)
{
if(ge->mesh_vertices.empty()) return;
if(ge->meshAttributes.method == MESH_TRANSFINITE) return;
bool forceOdd = false;
if((ge->meshAttributes.method != MESH_TRANSFINITE ||
CTX::instance()->mesh.flexibleTransfinite) &&
CTX::instance()->mesh.algoRecombine != 0) {
if(CTX::instance()->mesh.recombineAll) {
forceOdd = true;
}
}
if(!ge->getBeginVertex() || !ge->getEndVertex()) return;
MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
std::vector<std::pair<double, MVertex *> > lengths;
for(std::size_t i = 0; i < ge->mesh_vertices.size(); i++) {
MEdgeVertex *v = dynamic_cast<MEdgeVertex *>(ge->mesh_vertices[i]);
if(!v) {
Msg::Error("in 1D mesh filterPoints");
return;
}
double d = distance(v, v0);
double t;
v->getParameter(0, t);
if(i != 0) {
double t0;
if(v0->onWhat()->dim() == 0) {
// Vertex is begin point
t0 = ge->parFromPoint(SPoint3(v0->x(), v0->y(), v0->z()));
}
else
v0->getParameter(0, t0);
t = 0.5 * (t + t0);
}
double lc = F_LcB()(ge, t);
// double lc = v->getLc();
if(d < lc * .3) {
lengths.push_back(std::make_pair(lc / d, v));
}
else
v0 = v;
}
std::sort(lengths.begin(), lengths.end());
int last = lengths.size();
if(forceOdd) {
while(last % 2 != 0) last--;
}
/*
if(CTX::instance()->mesh.algoRecombine == 2){
if(last < 4) last = 0;
while (last %4 != 0)last--;
}
else{
while (last %2 != 0)last--;
}
}
*/
bool filteringObservesMinimumN =
(((int)ge->mesh_vertices.size() - last) >= nMinimumPoints);
if(filteringObservesMinimumN) {
for(int i = 0; i < last; i++) {
std::vector<MVertex *>::iterator it = std::find(
ge->mesh_vertices.begin(), ge->mesh_vertices.end(), lengths[i].second);
if(it != ge->mesh_vertices.end()) {
ge->mesh_vertices.erase(it);
}
delete lengths[i].second;
}
}
}
static void createPoints(GVertex *gv, GEdge *ge, BoundaryLayerField *blf,
std::vector<MVertex *> &v, const SVector3 &dir)
{
if(!ge->getBeginVertex() || !ge->getEndVertex()) return;
const double hwall = blf->hwall(gv->tag());
double L = hwall;
double LEdge = distance(ge->getBeginVertex()->mesh_vertices[0],
ge->getEndVertex()->mesh_vertices[0]);
while(1) {
if(L > blf->thickness || L > LEdge * .4) {
break;
}
SPoint3 p(gv->x() + dir.x() * L, gv->y() + dir.y() * L, 0.0);
v.push_back(new MEdgeVertex(p.x(), p.y(), p.z(), ge, ge->parFromPoint(p), 0,
blf->hfar));
int ith = v.size();
L += hwall * std::pow(blf->ratio, ith);
}
}
static void addBoundaryLayerPoints(GEdge *ge, double &t_begin, double &t_end,
std::vector<MVertex *> &_addBegin,
std::vector<MVertex *> &_addEnd)
{
_addBegin.clear();
_addEnd.clear();
// t_begin/t_end may change the left/right parameter of the interval
// _addBegin/_addEnd : additional points @ left/right
FieldManager *fields = ge->model()->getFields();
int n = fields->getNumBoundaryLayerFields();
if(n == 0) return;
if(!ge->getBeginVertex() || !ge->getEndVertex()) return;
// Check if edge is a BL edge
for(int i = 0; i < n; ++i) {
Field *bl_field = fields->get(fields->getBoundaryLayerField(i));
if(bl_field == NULL) continue;
BoundaryLayerField *blf = dynamic_cast<BoundaryLayerField *>(bl_field);
if(blf->isEdgeBL(ge->tag())) return;
}
SVector3 dir(ge->getEndVertex()->x() - ge->getBeginVertex()->x(),
ge->getEndVertex()->y() - ge->getBeginVertex()->y(),
ge->getEndVertex()->z() - ge->getBeginVertex()->z());
dir.normalize();
GVertex *gvb = ge->getBeginVertex();
GVertex *gve = ge->getEndVertex();
// Check if extremity nodes are BL node
for(int i = 0; i < n; ++i) {
Field *bl_field = fields->get(fields->getBoundaryLayerField(i));
if(bl_field == NULL) continue;
BoundaryLayerField *blf = dynamic_cast<BoundaryLayerField *>(bl_field);
blf->setupFor1d(ge->tag());
if(blf->isEndNode(gvb->tag())) {
if(ge->geomType() != GEntity::Line) {
Msg::Error("Boundary layer end point %d should lie on a straight line",
gvb->tag());
return;
}
if(_addBegin.size()) {
Msg::Error("Different boundary layers cannot touch each other");
return;
}
createPoints(gvb, ge, blf, _addBegin, dir);
if(!_addBegin.empty())
_addBegin[_addBegin.size() - 1]->getParameter(0, t_begin);
}
if(blf->isEndNode(gve->tag())) {
if(ge->geomType() != GEntity::Line) {
Msg::Error("Boundary layer end point %d should lie on a straight line",
gve->tag());
return;
}
if(_addEnd.size()) {
Msg::Error("Different boundary layers cannot touch each other");
return;
}
createPoints(gve, ge, blf, _addEnd, dir * -1.0);
if(!_addEnd.empty()) _addEnd[_addEnd.size() - 1]->getParameter(0, t_end);
}
}
}
void meshGEdge::operator()(GEdge *ge)
{
// debug stuff
if(CTX::instance()->debugSurface > 0){
std::vector<GFace *> f = ge->faces();
bool found = false;
for (size_t i=0;i<f.size(); i++){
if(f[i]->tag() == CTX::instance()->debugSurface) {
found = true;
}
}
if(!found) return;
}
ge->model()->setCurrentMeshEntity(ge);
if(ge->degenerate(1)) return;
if(ge->geomType() == GEntity::BoundaryLayerCurve) return;
if(ge->meshAttributes.method == MESH_NONE) return;
if(CTX::instance()->mesh.meshOnlyVisible && !ge->getVisibility()) return;
deMeshGEdge dem;
dem(ge);
if(MeshExtrudedCurve(ge)) return;
if(ge->getMeshMaster() != ge) {
GEdge *gef = dynamic_cast<GEdge *>(ge->getMeshMaster());
if(gef->meshStatistics.status == GEdge::PENDING) return;
Msg::Info("Meshing curve %d (%s) as a copy of %d", ge->tag(),
ge->getTypeString().c_str(), ge->getMeshMaster()->tag());
copyMesh(gef, ge, ge->masterOrientation);
ge->meshStatistics.status = GEdge::DONE;
return;
}
if(ge->model()->getNumEdges() > 100000) {
if(ge->tag() % 100000 == 1) {
Msg::Info("Meshing curve %d/%d (%s)", ge->tag(),
ge->model()->getNumEdges(), ge->getTypeString().c_str());
}
}
else if(ge->model()->getNumEdges() > 10000) {
if(ge->tag() % 10000 == 1) {
Msg::Info("Meshing curve %d/%d (%s)", ge->tag(),
ge->model()->getNumEdges(), ge->getTypeString().c_str());
}
}
else if(ge->model()->getNumEdges() > 1000) {
if(ge->tag() % 1000 == 1) {
Msg::Info("Meshing curve %d/%d (%s)", ge->tag(),
ge->model()->getNumEdges(), ge->getTypeString().c_str());
}
}
else {
Msg::Info("Meshing curve %d (%s)", ge->tag(), ge->getTypeString().c_str());
}
// compute bounds
Range<double> bounds = ge->parBounds(0);
double t_begin = bounds.low();
double t_end = bounds.high();
// if a BL is ending at one of the ends, then create specific points
std::vector<MVertex *> _addBegin, _addEnd;
addBoundaryLayerPoints(ge, t_begin, t_end, _addBegin, _addEnd);
// first compute the length of the curve by integrating one
double length;
std::vector<IntPoint> Points;
if(ge->geomType() == GEntity::Line && ge->getBeginVertex() &&
ge->getBeginVertex() == ge->getEndVertex() &&
// do not consider closed lines as degenerated
(ge->position(0.5) - ge->getBeginVertex()->xyz()).norm() <
CTX::instance()->geom.tolerance)
length = 0.0; // special case to avoid infinite loop in integration
else
length = Integration(ge, t_begin, t_end, F_One(), Points,
CTX::instance()->mesh.lcIntegrationPrecision *
CTX::instance()->lc);
ge->setLength(length);
Points.clear();
if(length < CTX::instance()->mesh.toleranceEdgeLength) {
ge->setTooSmall(true);
}
// Integrate detJ/lc du
double a;
int N;
int filterMinimumN = 1;
if(length == 0. && CTX::instance()->mesh.toleranceEdgeLength == 0.) {
Msg::Debug("Curve %d has a zero length", ge->tag());
a = 0.;
N = 1;
}
else if(ge->degenerate(0)) {
Msg::Debug("Curve %d is degenerated", ge->tag());
a = 0.;
N = 1;
}
else if(ge->meshAttributes.method == MESH_TRANSFINITE) {
a = Integration(ge, t_begin, t_end, F_Transfinite(), Points,
CTX::instance()->mesh.lcIntegrationPrecision);
N = ge->meshAttributes.nbPointsTransfinite;
if(CTX::instance()->mesh.flexibleTransfinite &&
CTX::instance()->mesh.lcFactor)
N /= CTX::instance()->mesh.lcFactor;
}
else {
if(CTX::instance()->mesh.algo2d == ALGO_2D_BAMG /* || blf*/) {
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);
}
// we should maybe provide an option to disable the smoothing
for(std::size_t i = 0; i < Points.size(); i++) {
IntPoint &pt = Points[i];
SVector3 der = ge->firstDer(pt.t);
pt.xp = der.norm();
}
if(CTX::instance()->mesh.algo2d != ALGO_2D_BAMG)
a = smoothPrimitive(ge, std::sqrt(CTX::instance()->mesh.smoothRatio),
Points);
filterMinimumN = ge->minimumMeshSegments() + 1;
N = std::max(filterMinimumN, (int)(a + 1.99));
}
// force odd number of points if blossom is used for recombination
if((ge->meshAttributes.method != MESH_TRANSFINITE ||
CTX::instance()->mesh.flexibleTransfinite) &&
CTX::instance()->mesh.algoRecombine != 0) {
std::vector<GFace *> const &faces = ge->faces();
if(CTX::instance()->mesh.recombineAll && faces.size()) {
if(N % 2 == 0) N++;
if(CTX::instance()->mesh.algoRecombine == 2) N = increaseN(N);
}
else {
for(std::vector<GFace *>::const_iterator it = faces.begin();
it != faces.end(); it++) {
if((*it)->meshAttributes.recombine) {
if(N % 2 == 0) N++;
if(CTX::instance()->mesh.algoRecombine == 2) N = increaseN(N);
break;
}
}
}
}
// printFandPrimitive(ge->tag(),Points);
// if the curve is periodic and if the begin vertex is identical to
// the end vertex and if this vertex has only one model curve
// adjacent to it, then the vertex is not connecting any other
// curve. So, the mesh vertex and its associated geom vertex are not
// necessary at the same location
// look if we are doing the STL triangulation
std::vector<MVertex *> &mesh_vertices = ge->mesh_vertices;
GPoint beg_p, end_p;
if(!ge->getBeginVertex() || !ge->getEndVertex()) {
Msg::Warning("Skipping curve with no begin or end point");
return;
}
else if(ge->getBeginVertex() == ge->getEndVertex() &&
ge->getBeginVertex()->edges().size() == 1) {
end_p = beg_p = ge->point(t_begin);
Msg::Debug("Meshing periodic closed curve");
}
else {
MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
MVertex *v1 = ge->getEndVertex()->mesh_vertices[0];
beg_p = GPoint(v0->x(), v0->y(), v0->z());
end_p = GPoint(v1->x(), v1->y(), v1->z());
}
// do not consider the first and the last vertex (those are not
// classified on this mesh edge)
if(N > 1) {
const double b = a / static_cast<double>(N - 1);
int count = 1, NUMP = 1;
IntPoint P1, P2;
mesh_vertices.resize(N - 2);
while(NUMP < N - 1) {
P1 = Points[count - 1];
P2 = Points[count];
const double d = (double)NUMP * b;
if((std::abs(P2.p) >= std::abs(d)) && (std::abs(P1.p) < std::abs(d))) {
double const dt = P2.t - P1.t;
double const dlc = P2.lc - P1.lc;
double const dp = P2.p - P1.p;
double const t = P1.t + dt / dp * (d - P1.p);
SVector3 der = ge->firstDer(t);
const double d = norm(der);
double lc = d / (P1.lc + dlc / dp * (d - P1.p));
GPoint V = ge->point(t);
mesh_vertices[NUMP - 1] =
new MEdgeVertex(V.x(), V.y(), V.z(), ge, t, 0, lc);
NUMP++;
}
else {
count++;
}
}
mesh_vertices.resize(NUMP - 1);
}
// Boundary Layer points are added
{
std::vector<MVertex *> vv;
vv.insert(vv.end(), _addBegin.begin(), _addBegin.end());
vv.insert(vv.end(), mesh_vertices.begin(), mesh_vertices.end());
for(std::size_t i = 0; i < _addEnd.size(); i++)
vv.push_back(_addEnd[_addEnd.size() - 1 - i]);
// vv.insert(vv.end(), _addEnd.rend(), _addEnd.rbegin());
mesh_vertices = vv;
}
if(CTX::instance()->mesh.algo2d != ALGO_2D_BAMG)
if(_addBegin.empty() && _addEnd.empty())
filterPoints(ge, filterMinimumN - 2);
std::vector<MLine *> &lines = ge->lines;
for(std::size_t i = 0; i <= mesh_vertices.size(); i++) {
MVertex *v0 =
(i == 0) ? ge->getBeginVertex()->mesh_vertices[0] : mesh_vertices[i - 1];
MVertex *v1 = (i == mesh_vertices.size()) ?
ge->getEndVertex()->mesh_vertices[0] :
mesh_vertices[i];
lines.push_back(new MLine(v0, v1));
}
if(ge->getBeginVertex() == ge->getEndVertex() &&
ge->getBeginVertex()->edges().size() == 1) {
MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
v0->x() = beg_p.x();
v0->y() = beg_p.y();
v0->z() = beg_p.z();
}
ge->meshStatistics.status = GEdge::DONE;
}
void orientMeshGEdge::operator()(GEdge *ge)
{
// apply user-specified mesh orientation constraints
if(ge->meshAttributes.reverseMesh)
for(std::size_t k = 0; k < ge->getNumMeshElements(); k++)
ge->getMeshElement(k)->reverse();
}