gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
gmshFace.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 
6 #include <stdlib.h>
7 #include "GModel.h"
8 #include "gmshFace.h"
9 #include "Geo.h"
10 #include "GeoInterpolation.h"
11 #include "Numeric.h"
12 #include "GmshMessage.h"
13 #include "Context.h"
14 #include "MTriangle.h"
15 #include "VertexArray.h"
16 
17 #if defined(HAVE_MESH)
18 #include "meshGFace.h"
19 #include "meshGEdge.h"
20 #include "Field.h"
21 #endif
22 
23 gmshFace::gmshFace(GModel *m, Surface *s) : GFace(m, s->Num)
24 {
25  resetNativePtr(s);
27 }
28 
29 bool gmshFace::degenerate(int dim) const
30 {
31  std::vector<GEdge *> const &eds = edges();
32  int numNonDegenerate = 0;
33  for(auto ge : eds) {
34  if(!ge->degenerate(0)) numNonDegenerate++;
35  }
36  return numNonDegenerate <= 1;
37 }
38 
40 {
41  _s = s;
42  l_edges.clear();
43  l_dirs.clear();
44  edgeLoops.clear();
45 
46  std::vector<GEdge *> eds;
47  std::vector<int> nums;
48  for(int i = 0; i < List_Nbr(_s->Generatrices); i++) {
49  Curve *c;
50  List_Read(_s->Generatrices, i, &c);
51  GEdge *e = model()->getEdgeByTag(abs(c->Num));
52  if(e) {
53  eds.push_back(e);
54  nums.push_back(c->Num);
55  }
56  else
57  Msg::Error("Unknown curve %d", c->Num);
58  }
59  for(int i = 0; i < List_Nbr(_s->GeneratricesByTag); i++) {
60  int j;
62  GEdge *e = model()->getEdgeByTag(abs(j));
63  if(e) {
64  eds.push_back(e);
65  nums.push_back(j);
66  }
67  else
68  Msg::Error("Unknown curve %d", j);
69  }
70 
71  std::vector<GEdge *> l_wire;
72  l_wire.reserve(eds.size());
73 
74  GVertex *first = nullptr;
75  for(std::size_t i = 0; i < eds.size(); i++) {
76  GEdge *e = eds[i];
77  int num = nums[i];
78  GVertex *start = (num > 0) ? e->getBeginVertex() : e->getEndVertex();
79  GVertex *next = (num > 0) ? e->getEndVertex() : e->getBeginVertex();
80  if(!first) first = start;
81  l_wire.push_back(e);
82  if(next == first) {
83  edgeLoops.push_back(GEdgeLoop(l_wire));
84  l_wire.clear();
85  first = nullptr;
86  }
87  l_edges.push_back(e);
88  e->addFace(this);
89  l_dirs.push_back((num > 0) ? 1 : -1);
90  if(List_Nbr(_s->Generatrices) == 2) {
92  std::max(e->meshAttributes.minimumMeshSegments, 2);
93  }
94  }
95 
96  // always compute and store the mean plane for plane surfaces (using
97  // the bounding vertices)
99 }
100 
102 {
103  if(!_s->geometry) return 1;
104  return _s->geometry->getMetricEigenvalue(pt);
105 }
106 
108 {
109  meshAttributes.recombine = _s->Recombine;
110  meshAttributes.recombineAngle = _s->RecombineAngle;
111  meshAttributes.method = _s->Method;
112  meshAttributes.extrude = _s->Extrude;
113  if(meshAttributes.method == MESH_TRANSFINITE) {
114  meshAttributes.transfiniteArrangement = _s->Recombine_Dir;
115  meshAttributes.transfiniteSmoothing = _s->TransfiniteSmoothing;
116  meshAttributes.corners.clear();
117  for(int i = 0; i < List_Nbr(_s->TrsfPoints); i++) {
118  Vertex *corn;
119  List_Read(_s->TrsfPoints, i, &corn);
120  GVertex *gv = model()->getVertexByTag(corn->Num);
121  if(gv)
122  meshAttributes.corners.push_back(gv);
123  else
124  Msg::Error("Unknown point %d in transfinite attributes", corn->Num);
125  }
126  }
127  meshAttributes.reverseMesh = _s->ReverseMesh;
128  meshAttributes.algorithm = _s->MeshAlgorithm;
129  meshAttributes.meshSizeFromBoundary = _s->MeshSizeFromBoundary;
130  meshAttributes.transfinite3 = false;
131 }
132 
133 Range<double> gmshFace::parBounds(int i) const { return Range<double>(0, 1); }
134 
135 SVector3 gmshFace::normal(const SPoint2 &param) const
136 {
137  if(_s->Typ != MSH_SURF_PLAN) {
138  Vertex vu = InterpolateSurface(_s, param[0], param[1], 1, 1);
139  Vertex vv = InterpolateSurface(_s, param[0], param[1], 1, 2);
140  Vertex n = vu % vv;
141  n.norme();
142  return SVector3(n.Pos.X, n.Pos.Y, n.Pos.Z);
143  }
144  else {
145  // We cannot use InterpolateSurface() for plane surfaces since it
146  // relies on the mean plane, which does not respect the
147  // orientation
148 
149  // FIXME: move this test at the end of the MeanPlane computation
150  // routine--and store the correct normal, damn it!
151 
152  double n[3] = {meanPlane.a, meanPlane.b, meanPlane.c};
153  norme(n);
154  GPoint pt = point(param.x(), param.y());
155  double v[3] = {pt.x(), pt.y(), pt.z()};
156  int NP = 10, tries = 0;
157  while(1) {
158  tries++;
159  double angle = 0.;
160  for(int i = 0; i < List_Nbr(_s->Generatrices); i++) {
161  Curve *c;
162  List_Read(_s->Generatrices, i, &c);
163  int N = (c->Typ == MSH_SEGM_LINE && List_Nbr(c->Control_Points) == 2) ?
164  1 :
165  NP;
166  for(int j = 0; j < N; j++) {
167  double u1 = (double)j / (double)N;
168  double u2 = (double)(j + 1) / (double)N;
169  Vertex p1 = InterpolateCurve(c, u1, 0);
170  Vertex p2 = InterpolateCurve(c, u2, 0);
171  double v1[3] = {p1.Pos.X, p1.Pos.Y, p1.Pos.Z};
172  double v2[3] = {p2.Pos.X, p2.Pos.Y, p2.Pos.Z};
173  angle += angle_plan(v, v1, v2, n);
174  }
175  }
176  if(fabs(angle) < 0.5) { // we're outside
177  NP *= 2;
178  Msg::Debug(
179  "Could not compute normal of surface %d - retrying with %d points",
180  tag(), NP);
181  if(tries > 10) {
182  Msg::Warning("Could not orient normal of surface %d", tag());
183  return SVector3(n[0], n[1], n[2]);
184  }
185  }
186  else if(angle > 0)
187  return SVector3(n[0], n[1], n[2]);
188  else
189  return SVector3(-n[0], -n[1], -n[2]);
190  }
191  }
192 }
193 
195 {
196  if(_s->Typ == MSH_SURF_PLAN && !_s->geometry) {
197  double x, y, z, VX[3], VY[3];
198  getMeanPlaneData(VX, VY, x, y, z);
199  return Pair<SVector3, SVector3>(SVector3(VX[0], VX[1], VX[2]),
200  SVector3(VY[0], VY[1], VY[2]));
201  }
202  else {
203  Vertex vu = InterpolateSurface(_s, param[0], param[1], 1, 1);
204  Vertex vv = InterpolateSurface(_s, param[0], param[1], 1, 2);
205  return Pair<SVector3, SVector3>(SVector3(vu.Pos.X, vu.Pos.Y, vu.Pos.Z),
206  SVector3(vv.Pos.X, vv.Pos.Y, vv.Pos.Z));
207  }
208 }
209 
210 void gmshFace::secondDer(const SPoint2 &param, SVector3 &dudu, SVector3 &dvdv,
211  SVector3 &dudv) const
212 {
213  if(_s->Typ == MSH_SURF_PLAN && !_s->geometry) {
214  dudu = SVector3(0., 0., 0.);
215  dvdv = SVector3(0., 0., 0.);
216  dudv = SVector3(0., 0., 0.);
217  }
218  else {
219  Vertex vuu = InterpolateSurface(_s, param[0], param[1], 2, 1);
220  Vertex vvv = InterpolateSurface(_s, param[0], param[1], 2, 2);
221  Vertex vuv = InterpolateSurface(_s, param[0], param[1], 2, 3);
222  dudu = SVector3(vuu.Pos.X, vuu.Pos.Y, vuu.Pos.Z);
223  dvdv = SVector3(vvv.Pos.X, vvv.Pos.Y, vvv.Pos.Z);
224  dudv = SVector3(vuv.Pos.X, vuv.Pos.Y, vuv.Pos.Z);
225  }
226 }
227 
228 GPoint gmshFace::point(double par1, double par2) const
229 {
230  double pp[2] = {par1, par2};
231  if(_s->Typ == MSH_SURF_PLAN && !_s->geometry) {
232  double x, y, z, VX[3], VY[3];
233  getMeanPlaneData(VX, VY, x, y, z);
234  return GPoint(x + VX[0] * par1 + VY[0] * par2,
235  y + VX[1] * par1 + VY[1] * par2,
236  z + VX[2] * par1 + VY[2] * par2, this, pp);
237  }
238  else {
239  Vertex v = InterpolateSurface(_s, par1, par2, 0, 0);
240  return GPoint(v.Pos.X, v.Pos.Y, v.Pos.Z, this, pp);
241  }
242 }
243 
245  const double initialGuess[2]) const
246 {
247 #if defined(HAVE_ALGLIB)
248  return GFace::closestPoint(qp, initialGuess);
249 #endif
250  if(_s->Typ == MSH_SURF_PLAN && !_s->geometry) {
251  double XP = qp.x();
252  double YP = qp.y();
253  double ZP = qp.z();
254  double VX[3], VY[3], x, y, z;
255  getMeanPlaneData(VX, VY, x, y, z);
256  double M[3][2] = {{VX[0], VY[0]}, {VX[1], VY[1]}, {VX[2], VY[2]}};
257  double MN[2][2];
258  double B[3] = {XP - x, YP - y, ZP - z};
259  double BN[2], UV[2];
260  for(int i = 0; i < 2; i++) {
261  BN[i] = 0;
262  for(int k = 0; k < 3; k++) { BN[i] += B[k] * M[k][i]; }
263  }
264  for(int i = 0; i < 2; i++) {
265  for(int j = 0; j < 2; j++) {
266  MN[i][j] = 0;
267  for(int k = 0; k < 3; k++) { MN[i][j] += M[k][i] * M[k][j]; }
268  }
269  }
270  sys2x2(MN, BN, UV);
271  return GPoint(XP, YP, ZP, this, UV);
272  }
273 
274  Vertex v;
275  v.Pos.X = qp.x();
276  v.Pos.Y = qp.y();
277  v.Pos.Z = qp.z();
278  double u[2] = {initialGuess[0], initialGuess[1]};
279  bool result = ProjectPointOnSurface(_s, v, u);
280  if(!result) return GPoint(-1.e22, -1.e22, -1.e22, nullptr, u);
281  return GPoint(v.Pos.X, v.Pos.Y, v.Pos.Z, this, u);
282 }
283 
284 SPoint2 gmshFace::parFromPoint(const SPoint3 &qp, bool onSurface,
285  bool convTestXYZ) const
286 {
287  if(_s->Typ == MSH_SURF_PLAN) {
288  double x, y, z, VX[3], VY[3];
289  getMeanPlaneData(VX, VY, x, y, z);
290  double const vec[3] = {qp.x() - x, qp.y() - y, qp.z() - z};
291  double const u = prosca(vec, VX);
292  double const v = prosca(vec, VY);
293  return SPoint2(u, v);
294  }
295  else {
296  return GFace::parFromPoint(qp, onSurface, convTestXYZ);
297  }
298 }
299 
301 {
302  switch(_s->Typ) {
303  case MSH_SURF_PLAN:
304  if(_s->geometry)
305  return ParametricSurface;
306  else
307  return Plane;
308  case MSH_SURF_REGL:
309  case MSH_SURF_TRIC: return RuledSurface;
310  case MSH_SURF_DISCRETE: return DiscreteSurface;
312  default: return Unknown;
313  }
314 }
315 
317 {
318  return geomType() != BoundaryLayerSurface;
319 }
320 
321 bool gmshFace::containsPoint(const SPoint3 &pt) const
322 {
323  if(_s->Typ == MSH_SURF_PLAN) {
324  // OK to use the normal from the mean plane here: we compensate for the
325  // (possibly wrong) orientation at the end
326  double n[3] = {meanPlane.a, meanPlane.b, meanPlane.c};
327  norme(n);
328  double angle = 0.;
329  double v[3] = {pt.x(), pt.y(), pt.z()};
330  for(int i = 0; i < List_Nbr(_s->Generatrices); i++) {
331  Curve *c;
332  List_Read(_s->Generatrices, i, &c);
333  int N = (c->Typ == MSH_SEGM_LINE) ? 1 : 10;
334  for(int j = 0; j < N; j++) {
335  double u1 = (double)j / (double)N;
336  double u2 = (double)(j + 1) / (double)N;
337  Vertex p1 = InterpolateCurve(c, u1, 0);
338  Vertex p2 = InterpolateCurve(c, u2, 0);
339  double v1[3] = {p1.Pos.X, p1.Pos.Y, p1.Pos.Z};
340  double v2[3] = {p2.Pos.X, p2.Pos.Y, p2.Pos.Z};
341  angle += angle_plan(v, v1, v2, n);
342  }
343  }
344  // we're inside if angle equals 2 * pi
345  if(fabs(angle) > 2 * M_PI - 0.5 && fabs(angle) < 2 * M_PI + 0.5)
346  return true;
347  return false;
348  }
349 
350  return false;
351 }
gmshFace::parBounds
virtual Range< double > parBounds(int i) const
Definition: gmshFace.cpp:133
Surface::TrsfPoints
List_T * TrsfPoints
Definition: Geo.h:122
MTriangle.h
Geo.h
GeoInterpolation.h
GEntity::BoundaryLayerSurface
@ BoundaryLayerSurface
Definition: GEntity.h:116
GFace::edges
virtual std::vector< GEdge * > const & edges() const
Definition: GFace.h:121
Field.h
GPoint::y
double y() const
Definition: GPoint.h:22
mean_plane::c
double c
Definition: Numeric.h:35
gmshFace::geomType
virtual GEntity::GeomType geomType() const
Definition: gmshFace.cpp:300
Curve
Definition: Geo.h:74
gmshFace::_s
Surface * _s
Definition: gmshFace.h:15
Surface::Typ
int Typ
Definition: Geo.h:114
GFace
Definition: GFace.h:33
GEntity::ParametricSurface
@ ParametricSurface
Definition: GEntity.h:112
GEntity::model
GModel * model() const
Definition: GEntity.h:277
MSH_SURF_BND_LAYER
#define MSH_SURF_BND_LAYER
Definition: GeoDefines.h:36
angle
double angle(const SVector3 &a, const SVector3 &b)
Definition: SVector3.h:157
SPoint2
Definition: SPoint2.h:12
MSH_SURF_PLAN
#define MSH_SURF_PLAN
Definition: GeoDefines.h:33
Msg::Debug
static void Debug(const char *fmt,...)
Definition: GmshMessage.cpp:752
c
static double c(int i, int j, fullMatrix< double > &CA, const std::vector< SPoint3 > &P, const std::vector< SPoint3 > &Q)
Definition: discreteFrechetDistance.cpp:15
MESH_TRANSFINITE
#define MESH_TRANSFINITE
Definition: GmshDefines.h:259
Msg::Warning
static void Warning(const char *fmt,...)
Definition: GmshMessage.cpp:543
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
VertexArray.h
GEntity::Unknown
@ Unknown
Definition: GEntity.h:89
SPoint3
Definition: SPoint3.h:14
List_Nbr
int List_Nbr(List_T *liste)
Definition: ListUtils.cpp:106
Vertex
Definition: Geo.h:29
GFace::meshAttributes
struct GFace::@18 meshAttributes
SVector3
Definition: SVector3.h:16
gmshSurface::getMetricEigenvalue
virtual double getMetricEigenvalue(const SPoint2 &)
Definition: gmshSurface.cpp:32
Surface::MeshSizeFromBoundary
int MeshSizeFromBoundary
Definition: Geo.h:130
meshGEdge.h
mean_plane::b
double b
Definition: Numeric.h:35
GModel::getEdgeByTag
GEdge * getEdgeByTag(int n) const
Definition: GModel.cpp:336
GmshMessage.h
gmshFace::resetNativePtr
void resetNativePtr(Surface *s)
Definition: gmshFace.cpp:39
GPoint
Definition: GPoint.h:13
GFace::closestPoint
virtual GPoint closestPoint(const SPoint3 &queryPoint, const double initialGuess[2]) const
Definition: GFace.cpp:1323
GFace::l_dirs
std::vector< int > l_dirs
Definition: GFace.h:38
InterpolateCurve
Vertex InterpolateCurve(Curve *c, double u, int const derivee)
Definition: GeoInterpolation.cpp:441
GEdge::addFace
void addFace(GFace *f)
Definition: GEdge.cpp:200
GEntity::Plane
@ Plane
Definition: GEntity.h:105
Vertex::Num
int Num
Definition: Geo.h:31
GPoint::z
double z() const
Definition: GPoint.h:23
ProjectPointOnSurface
bool ProjectPointOnSurface(Surface *s, Vertex &p, double uv[2])
Definition: Geo.cpp:3340
gmshFace::firstDer
virtual Pair< SVector3, SVector3 > firstDer(const SPoint2 &param) const
Definition: gmshFace.cpp:194
SPoint3::x
double x(void) const
Definition: SPoint3.h:125
GEdgeLoop
Definition: GEdgeLoop.h:36
Range
Definition: Range.h:10
GFace::edgeLoops
std::vector< GEdgeLoop > edgeLoops
Definition: GFace.h:47
gmshFace.h
Surface::Extrude
ExtrudeParams * Extrude
Definition: Geo.h:124
angle_plan
double angle_plan(double v[3], double p1[3], double p2[3], double n[3])
Definition: Numeric.cpp:267
GVertex
Definition: GVertex.h:23
Vertex::norme
void norme()
Definition: Geo.h:51
gmshFace::resetMeshAttributes
virtual void resetMeshAttributes()
Definition: gmshFace.cpp:107
Numeric.h
gmshFace::getMetricEigenvalue
virtual double getMetricEigenvalue(const SPoint2 &)
Definition: gmshFace.cpp:101
GModel
Definition: GModel.h:44
GFace::getMeanPlaneData
void getMeanPlaneData(double VX[3], double VY[3], double &x, double &y, double &z) const
Definition: GFace.cpp:902
GEdge::getBeginVertex
virtual GVertex * getBeginVertex() const
Definition: GEdge.h:63
InterpolateSurface
Vertex InterpolateSurface(Surface *s, double u, double v, int derivee, int u_v)
Definition: GeoInterpolation.cpp:958
GFace::parFromPoint
virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true, bool convTestXYZ=false) const
Definition: GFace.cpp:1254
Surface
Definition: Geo.h:111
MSH_SURF_DISCRETE
#define MSH_SURF_DISCRETE
Definition: GeoDefines.h:38
mean_plane::a
double a
Definition: Numeric.h:35
gmshFace::closestPoint
virtual GPoint closestPoint(const SPoint3 &queryPoint, const double initialGuess[2]) const
Definition: gmshFace.cpp:244
GEdge::meshAttributes
struct GEdge::@16 meshAttributes
SPoint3::y
double y(void) const
Definition: SPoint3.h:127
gmshFace::haveParametrization
virtual bool haveParametrization()
Definition: gmshFace.cpp:316
prosca
double prosca(double const a[3], double const b[3])
Definition: Numeric.h:112
MSH_SURF_REGL
#define MSH_SURF_REGL
Definition: GeoDefines.h:34
GEntity::tag
int tag() const
Definition: GEntity.h:280
GEntity::GeomType
GeomType
Definition: GEntity.h:88
GEntity::DiscreteSurface
@ DiscreteSurface
Definition: GEntity.h:117
Surface::GeneratricesByTag
List_T * GeneratricesByTag
Definition: Geo.h:121
MSH_SURF_TRIC
#define MSH_SURF_TRIC
Definition: GeoDefines.h:35
gmshFace::secondDer
virtual void secondDer(const SPoint2 &, SVector3 &, SVector3 &, SVector3 &) const
Definition: gmshFace.cpp:210
Surface::MeshAlgorithm
int MeshAlgorithm
Definition: Geo.h:130
gmshFace::parFromPoint
virtual SPoint2 parFromPoint(const SPoint3 &, bool onSurface=true, bool convTestXYZ=false) const
Definition: gmshFace.cpp:284
Surface::TransfiniteSmoothing
int TransfiniteSmoothing
Definition: Geo.h:119
meshGFace.h
SPoint2::x
double x(void) const
Definition: SPoint2.h:86
Coord::X
double X
Definition: Geo.h:26
Surface::RecombineAngle
double RecombineAngle
Definition: Geo.h:118
Surface::Recombine
int Recombine
Definition: Geo.h:116
Context.h
sys2x2
int sys2x2(double mat[2][2], double b[2], double res[2])
Definition: Numeric.cpp:101
z
const double z
Definition: GaussQuadratureQuad.cpp:56
Pair
Definition: Pair.h:10
Coord::Z
double Z
Definition: Geo.h:26
GEdge
Definition: GEdge.h:26
GFace::meanPlane
mean_plane meanPlane
Definition: GFace.h:40
gmshFace::point
virtual GPoint point(double par1, double par2) const
Definition: gmshFace.cpp:228
gmshFace::containsPoint
virtual bool containsPoint(const SPoint3 &pt) const
Definition: gmshFace.cpp:321
Coord::Y
double Y
Definition: Geo.h:26
SPoint3::z
double z(void) const
Definition: SPoint3.h:129
Vertex::Pos
Coord Pos
Definition: Geo.h:34
Surface::Generatrices
List_T * Generatrices
Definition: Geo.h:120
GModel.h
Surface::geometry
gmshSurface * geometry
Definition: Geo.h:129
Surface::Recombine_Dir
int Recombine_Dir
Definition: Geo.h:117
GEdge::getEndVertex
virtual GVertex * getEndVertex() const
Definition: GEdge.h:64
gmshFace::normal
virtual SVector3 normal(const SPoint2 &param) const
Definition: gmshFace.cpp:135
Surface::ReverseMesh
int ReverseMesh
Definition: Geo.h:130
gmshFace::degenerate
bool degenerate(int dim) const
Definition: gmshFace.cpp:29
norme
double norme(double a[3])
Definition: Numeric.h:123
GFace::l_edges
std::vector< GEdge * > l_edges
Definition: GFace.h:37
GModel::getVertexByTag
GVertex * getVertexByTag(int n) const
Definition: GModel.cpp:346
gmshFace::gmshFace
gmshFace(GModel *m, Surface *s)
Definition: gmshFace.cpp:23
GPoint::x
double x() const
Definition: GPoint.h:21
GEdge::minimumMeshSegments
int minimumMeshSegments
Definition: GEdge.h:255
List_Read
void List_Read(List_T *liste, int index, void *data)
Definition: ListUtils.cpp:111
MSH_SEGM_LINE
#define MSH_SEGM_LINE
Definition: GeoDefines.h:20
GFace::computeMeanPlane
void computeMeanPlane()
Definition: GFace.cpp:645
GEntity::RuledSurface
@ RuledSurface
Definition: GEntity.h:111
Surface::Method
int Method
Definition: Geo.h:115
SPoint2::y
double y(void) const
Definition: SPoint2.h:88