1 // MeshOptimizer - Copyright (C) 2013-2019 UCLouvain-ULiege
2 //
3 // Permission is hereby granted, free of charge, to any person
4 // obtaining a copy of this software and associated documentation
5 // files (the "Software"), to deal in the Software without
6 // restriction, including without limitation the rights to use, copy,
7 // modify, merge, publish, distribute, and/or sell copies of the
8 // Software, and to permit persons to whom the Software is furnished
9 // to do so, provided that the above copyright notice(s) and this
10 // permission notice appear in all copies of the Software and that
11 // both the above copyright notice(s) and this permission notice
12 // appear in supporting documentation.
13 //
14 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
16 // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17 // NONINFRINGEMENT OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE
18 // COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR
19 // ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY
20 // DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
21 // WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
22 // ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
23 // OF THIS SOFTWARE.
24 //
25 // Contributor(s): Thomas Toulorge, Jonathan Lambrechts
26
27 #include <iostream>
28 #include <algorithm>
29 #include "GmshMessage.h"
30 #include "GFace.h"
31 #include "GEdge.h"
32 #include "MVertex.h"
33 #include "MLine.h"
34 #include "VertexCoord.h"
35
getUvw(MVertex * vert) const36 SPoint3 VertexCoordParent::getUvw(MVertex *vert) const
37 {
38 GEntity *ge = vert->onWhat();
39 switch(ge->dim()) {
40 case 1: {
41 SPoint3 p(0., 0., 0.);
42 reparamMeshVertexOnEdge(
43 vert, static_cast<GEdge *>(ge),
44 p[0]); // Overkill if vert. well classified and parametrized
45 return p;
46 break;
47 }
48 case 2: {
49 SPoint2 p;
50 reparamMeshVertexOnFace(
51 vert, static_cast<GFace *>(ge),
52 p); // Overkill if vert. well classified and parametrized
53 return SPoint3(p[0], p[1], 0.);
54 break;
55 }
56 }
57 return SPoint3(0., 0., 0.);
58 }
59
uvw2Xyz(const SPoint3 & uvw) const60 SPoint3 VertexCoordParent::uvw2Xyz(const SPoint3 &uvw) const
61 {
62 GEntity *ge = _vert->onWhat();
63 GPoint gp = (ge->dim() == 1) ?
64 static_cast<GEdge *>(ge)->point(uvw[0]) :
65 static_cast<GFace *>(ge)->point(uvw[0], uvw[1]);
66 return SPoint3(gp.x(), gp.y(), gp.z());
67 }
68
gXyz2gUvw(const SPoint3 & uvw,const SPoint3 & gXyz,SPoint3 & gUvw) const69 void VertexCoordParent::gXyz2gUvw(const SPoint3 &uvw, const SPoint3 &gXyz,
70 SPoint3 &gUvw) const
71 {
72 GEntity *ge = _vert->onWhat();
73
74 if(ge->dim() == 1) {
75 SVector3 der = static_cast<GEdge *>(ge)->firstDer(uvw[0]);
76 gUvw[0] = gXyz.x() * der.x() + gXyz.y() * der.y() + gXyz.z() * der.z();
77 }
78 else {
79 Pair<SVector3, SVector3> der =
80 static_cast<GFace *>(ge)->firstDer(SPoint2(uvw[0], uvw[1]));
81 gUvw[0] = gXyz.x() * der.first().x() + gXyz.y() * der.first().y() +
82 gXyz.z() * der.first().z();
83 gUvw[1] = gXyz.x() * der.second().x() + gXyz.y() * der.second().y() +
84 gXyz.z() * der.second().z();
85 }
86 }
87
gXyz2gUvw(const SPoint3 & uvw,const std::vector<SPoint3> & gXyz,std::vector<SPoint3> & gUvw) const88 void VertexCoordParent::gXyz2gUvw(const SPoint3 &uvw,
89 const std::vector<SPoint3> &gXyz,
90 std::vector<SPoint3> &gUvw) const
91 {
92 GEntity *ge = _vert->onWhat();
93
94 if(ge->dim() == 1) {
95 SVector3 der = static_cast<GEdge *>(ge)->firstDer(uvw[0]);
96 auto itUvw = gUvw.begin();
97 for(auto itXyz = gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
98 (*itUvw)[0] =
99 itXyz->x() * der.x() + itXyz->y() * der.y() + itXyz->z() * der.z();
100 itUvw++;
101 }
102 }
103 else {
104 Pair<SVector3, SVector3> der =
105 static_cast<GFace *>(ge)->firstDer(SPoint2(uvw[0], uvw[1]));
106 auto itUvw = gUvw.begin();
107 for(auto itXyz = gXyz.begin(); itXyz != gXyz.end(); itXyz++) {
108 (*itUvw)[0] = itXyz->x() * der.first().x() +
109 itXyz->y() * der.first().y() + itXyz->z() * der.first().z();
110 (*itUvw)[1] = itXyz->x() * der.second().x() +
111 itXyz->y() * der.second().y() +
112 itXyz->z() * der.second().z();
113 itUvw++;
114 }
115 }
116 }
117
118 namespace {
119
getLineElTangent(MElement * el,int iNode)120 SVector3 getLineElTangent(MElement *el, int iNode)
121 {
122 double gsf[1256][3], u, v, w;
123 el->getNode(iNode, u, v, w);
124 // el->getGradShapeFunctions(u,v,w,gsf);
125 el->getGradShapeFunctions(u, v, w, gsf, 1);
126
127 SVector3 dxyzdu(0.);
128 // int nSF = el->getNumShapeFunctions()();
129 int nSF = el->getNumPrimaryVertices();
130 for(int j = 0; j < nSF; j++) {
131 const SPoint3 p = el->getVertex(j)->point();
132 dxyzdu(0) += gsf[j][0] * p.x();
133 dxyzdu(1) += gsf[j][0] * p.y();
134 dxyzdu(2) += gsf[j][0] * p.z();
135 }
136 dxyzdu.normalize();
137
138 return dxyzdu;
139 }
140
getSurfElNormal(MElement * el,int iNode)141 SVector3 getSurfElNormal(MElement *el, int iNode)
142 {
143 double gsf[1256][3], u, v, w;
144 el->getNode(iNode, u, v, w);
145 // el->getGradShapeFunctions(u,v,w,gsf);
146 el->getGradShapeFunctions(u, v, w, gsf, 1);
147
148 SVector3 dxyzdu(0.), dxyzdv(0.);
149 // int nSF = el->getNumShapeFunctions()();
150 int nSF = el->getNumPrimaryVertices();
151 for(int j = 0; j < nSF; j++) {
152 const SPoint3 p = el->getVertex(j)->point();
153 dxyzdu(0) += gsf[j][0] * p.x();
154 dxyzdu(1) += gsf[j][0] * p.y();
155 dxyzdu(2) += gsf[j][0] * p.z();
156 dxyzdv(0) += gsf[j][1] * p.x();
157 dxyzdv(1) += gsf[j][1] * p.y();
158 dxyzdv(2) += gsf[j][1] * p.z();
159 }
160
161 SVector3 normal = crossprod(dxyzdu, dxyzdv);
162 normal.normalize();
163 return normal;
164 }
165
166 } // namespace
167
VertexCoordLocalLine(MVertex * v)168 VertexCoordLocalLine::VertexCoordLocalLine(MVertex *v)
169 : dir(0.), x0(v->x()), y0(v->y()), z0(v->z())
170 {
171 GEntity *ge = v->onWhat();
172 const unsigned nEl = ge->getNumMeshElements();
173
174 for(unsigned iEl = 0; iEl < nEl; iEl++) {
175 MElement *el = ge->getMeshElement(iEl);
176 std::vector<MVertex *> lVerts;
177 el->getVertices(lVerts);
178 auto itV = std::find(lVerts.begin(), lVerts.end(), v);
179 if(itV != lVerts.end()) {
180 const int iNode = std::distance(lVerts.begin(), itV);
181 dir += getLineElTangent(el, iNode);
182 }
183 }
184 dir.normalize();
185 }
186
VertexCoordLocalSurf(MVertex * v)187 VertexCoordLocalSurf::VertexCoordLocalSurf(MVertex *v)
188 : x0(v->x()), y0(v->y()), z0(v->z())
189 {
190 GEntity *ge = v->onWhat();
191 const unsigned nEl = ge->getNumMeshElements();
192
193 SVector3 n(0.);
194 for(unsigned iEl = 0; iEl < nEl; iEl++) {
195 MElement *el = ge->getMeshElement(iEl);
196 std::vector<MVertex *> lVerts;
197 el->getVertices(lVerts);
198 auto itV = std::find(lVerts.begin(), lVerts.end(), v);
199 if(itV != lVerts.end()) {
200 const int iNode = std::distance(lVerts.begin(), itV);
201 n += getSurfElNormal(el, iNode);
202 }
203 }
204 n.normalize();
205
206 if(fabs(fabs(dot(n, SVector3(1., 0., 0.))) - 1.) <
207 1.e-10) { // If normal is x-axis, take y- and z- axis as dir.
208 dir0 = SVector3(0., 1., 0.);
209 dir1 = SVector3(0., 0., 1.);
210 }
211 else {
212 dir0 =
213 SVector3(1. - n.x() * n.x(), -n.x() * n.y(),
214 -n.x() * n.z()); // 1st dir. = (normalized) proj. of e_x in plane
215 dir0.normalize();
216 dir1 = crossprod(dir0, n);
217 }
218 }
219