Skip to content

Commit 250cbca

Browse files
Nightly Botcwsmith
Nightly Bot
authored and
cwsmith
committed
Merging develop into master
2 parents 2c4c41f + f3ce409 commit 250cbca

File tree

5 files changed

+445
-6
lines changed

5 files changed

+445
-6
lines changed

Diff for: ‎ma/CMakeLists.txt

+1
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@ set(SOURCES
2020
maSize.cc
2121
maOperator.cc
2222
maCollapse.cc
23+
maRegionCollapse.cc
2324
maMatchedCollapse.cc
2425
maLayerCollapse.cc
2526
maMatch.cc

Diff for: ‎ma/maRegionCollapse.cc

+385
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,385 @@
1+
/******************************************************************************
2+
3+
Copyright 2013 Scientific Computation Research Center,
4+
Rensselaer Polytechnic Institute. All rights reserved.
5+
6+
The LICENSE file included with this distribution describes the terms
7+
of the SCOREC Non-Commercial License this program is distributed under.
8+
9+
*******************************************************************************/
10+
#include "maRegionCollapse.h"
11+
#include "maAdapt.h"
12+
#include "maShape.h"
13+
#include <apfCavityOp.h>
14+
#include <pcu_util.h>
15+
16+
namespace ma {
17+
18+
void RegionCollapse::Init(Adapt* a, double fa)
19+
{
20+
adapt = a;
21+
region = 0;
22+
reclassifyEdge = 0;
23+
numBdryFaces = 0;
24+
faces[0] = faces[1] = faces[2] = faces[3] = 0;
25+
flatAngle = fa;
26+
}
27+
28+
bool RegionCollapse::requestLocality(apf::CavityOp* o)
29+
{
30+
/* get vertices again since this is sometimes used
31+
before setVerts can work */
32+
Entity* v[4];
33+
Mesh* m = adapt->mesh;
34+
m->getDownward(region,0,v);
35+
return o->requestLocality(v,4);
36+
}
37+
38+
bool RegionCollapse::setRegion(Entity* r)
39+
{
40+
if (getFlag(adapt,r,DONT_COLLAPSE))
41+
return false;
42+
region = r;
43+
return true;
44+
}
45+
46+
static int relativeDirection(Vector n1, Vector n2, double flatAngle)
47+
{
48+
double cosFlatAngleSquare;
49+
double dot,ll1,ll2;
50+
51+
cosFlatAngleSquare = cos(flatAngle*0.017453293); // PI/180
52+
cosFlatAngleSquare = cosFlatAngleSquare * cosFlatAngleSquare;
53+
54+
dot = n1 * n2;
55+
56+
ll1 = n1 * n1;
57+
ll2 = n2 * n2;
58+
59+
if(dot*dot/(ll1*ll2) < cosFlatAngleSquare) return 0;
60+
if(dot > 0) return 1;
61+
return -1;
62+
}
63+
64+
static Entity* getTetVertOppositeFace(Mesh* m, Entity* r, Entity* f)
65+
{
66+
Entity* rfs[4];
67+
m->getDownward(r, 2, rfs);
68+
int index = findIn(rfs, 4, f);
69+
PCU_ALWAYS_ASSERT(index > -1);
70+
71+
Entity* rvs[4];
72+
Entity* fvs[3];
73+
74+
m->getDownward(r, 0, rvs);
75+
m->getDownward(f, 0, fvs);
76+
77+
for (int i = 0; i < 4; i++) {
78+
int index = findIn(fvs, 3, rvs[i]);
79+
if (index == -1)
80+
return rvs[i];
81+
}
82+
return 0;
83+
}
84+
85+
Vector faceNormal(Mesh* mesh, Entity* face, Entity* rgn)
86+
{
87+
Vector xyz[4];
88+
Vector v01, v02, v03;
89+
Entity* fvs[3];
90+
mesh->getDownward(face, 0, fvs);
91+
for (int i = 0; i < 3; i++) {
92+
xyz[i] = getPosition(mesh, fvs[i]);
93+
}
94+
95+
96+
Entity* vertex = getTetVertOppositeFace(mesh, rgn, face);
97+
xyz[3] = getPosition(mesh, vertex);
98+
99+
v01 = xyz[1] - xyz[0];
100+
v02 = xyz[2] - xyz[0];
101+
v03 = xyz[3] - xyz[0];
102+
103+
Vector normal = cross(v01, v02);
104+
if( (normal * v03) > 0 )
105+
normal = normal * (-1);
106+
107+
return normal;
108+
}
109+
110+
bool RegionCollapse::checkGeom()
111+
{
112+
Vector normal1, normal2;
113+
Mesh* mesh = adapt->mesh;
114+
switch (numBdryFaces) {
115+
case 1:
116+
normal1 = faceNormal(mesh, faces[0], region);
117+
for (int i = 1; i < 4; i++) {
118+
normal2 = faceNormal(mesh, faces[i], region);
119+
if (relativeDirection(normal1, normal2, flatAngle) != -1)
120+
return false;
121+
}
122+
break;
123+
case 2:
124+
normal1 = faceNormal(mesh, faces[0], region);
125+
normal2 = faceNormal(mesh, faces[1], region);
126+
if (relativeDirection(normal1, normal2, flatAngle) != 1)
127+
return false;
128+
normal1 = faceNormal(mesh, faces[2], region);
129+
normal2 = faceNormal(mesh, faces[3], region);
130+
if (relativeDirection(normal1, normal2, flatAngle) != 1)
131+
return false;
132+
break;
133+
case 3:
134+
normal1 = faceNormal(mesh, faces[3], region);
135+
for (int i = 0; i < 3; i++) {
136+
normal2 = faceNormal(mesh, faces[i], region);
137+
if (!relativeDirection(normal1, normal2, flatAngle))
138+
return false;
139+
}
140+
break;
141+
default:
142+
return false;
143+
}
144+
return true;
145+
}
146+
147+
148+
bool RegionCollapse::checkTopo()
149+
{
150+
151+
Mesh* mesh = adapt->mesh;
152+
153+
Model* modelFace = 0;
154+
Model* uniqueModelFace = 0;
155+
numBdryFaces = 0;
156+
157+
Entity* fs[4];
158+
mesh->getDownward(region, 2, fs);
159+
int m = 0;
160+
for (int i = 0; i < 4; i++) {
161+
Entity* face = fs[i];
162+
163+
int faceModelType = mesh->getModelType(mesh->toModel(face));
164+
if (faceModelType != 2) {
165+
faces[3-m] = face;
166+
++m;
167+
continue;
168+
}
169+
170+
if (mesh->countUpward(face) == 2) {
171+
faces[3-m] = face;
172+
++m;
173+
continue;
174+
}
175+
176+
faces[numBdryFaces] = face;
177+
++numBdryFaces;
178+
179+
modelFace = mesh->toModel(face);
180+
if (uniqueModelFace) {
181+
if (uniqueModelFace != modelFace)
182+
return false;
183+
}
184+
else
185+
uniqueModelFace = modelFace;
186+
}
187+
188+
Entity* oppositeEdge = 0;
189+
switch (numBdryFaces) {
190+
case 0:
191+
return false;
192+
case 1:
193+
{
194+
Entity* vert = getTetVertOppositeFace(mesh, region, faces[0]);
195+
if (mesh->getModelType(mesh->toModel(vert)) != 3)
196+
return false;
197+
break;
198+
}
199+
case 2:
200+
{
201+
Entity* f2es[3];
202+
Entity* f3es[3];
203+
mesh->getDownward(faces[2], 1, f2es);
204+
mesh->getDownward(faces[3], 1, f3es);
205+
int j;
206+
for (j = 0; j < 3; j++) {
207+
int index = findIn(f3es, 3, f2es[j]);
208+
if (index > -1)
209+
break;
210+
}
211+
reclassifyEdge = f2es[j]; // this is the edge being reclassified
212+
213+
Entity* f0es[3];
214+
Entity* f1es[3];
215+
mesh->getDownward(faces[0], 1, f0es);
216+
mesh->getDownward(faces[1], 1, f1es);
217+
for (j = 0; j < 3; j++) {
218+
int index = findIn(f1es, 3, f0es[j]);
219+
if (index > -1)
220+
break;
221+
}
222+
oppositeEdge = f0es[j]; // this is the edge being deleted
223+
224+
if (mesh->getModelType(mesh->toModel(reclassifyEdge)) != 3)
225+
return false;
226+
if (mesh->getModelType(mesh->toModel(oppositeEdge)) == 1)
227+
return false;
228+
break;
229+
}
230+
case 3:
231+
{
232+
Entity* f3es[3];
233+
for (int j = 0; j < 2; j++) {
234+
Entity* fes[3];
235+
mesh->getDownward(faces[j], 1, fes);
236+
for (int k = 0; k < 3; k++) {
237+
int index = findIn(f3es, 3, fes[k]);
238+
if (index > -1) continue;
239+
if (mesh->getModelType(mesh->toModel(fes[k])) != 2)
240+
return false;
241+
}
242+
}
243+
Entity* vert = getTetVertOppositeFace(mesh, region, faces[3]);
244+
if (mesh->getModelType(mesh->toModel(vert)) != 2)
245+
return false;
246+
if (mesh->countUpward(vert) != 3)
247+
return false;
248+
break;
249+
}
250+
default:
251+
return false;
252+
}
253+
254+
return true;
255+
}
256+
257+
static void processNewBdryVert(Mesh* m, Entity* v)
258+
{
259+
Model* c = m->toModel(v);
260+
int modelType = m->getModelType(c);
261+
PCU_ALWAYS_ASSERT(modelType == 1 || modelType == 2);
262+
263+
Vector coords = getPosition(m, v);
264+
Vector newCoords;
265+
Vector newParams;
266+
m->getClosestPoint(c, coords, newCoords, newParams);
267+
m->setParam(v, newParams);
268+
// TODO: this must be done via snapping
269+
/* m->setPoint(v, 0, newCoords); */
270+
}
271+
272+
void RegionCollapse::apply()
273+
{
274+
Mesh* mesh = adapt->mesh;
275+
Model* modelFace = mesh->toModel(faces[0]);
276+
277+
switch (numBdryFaces) {
278+
case 1:
279+
{
280+
// reclassify faces
281+
for (int i = 1; i < 4; i++) {
282+
mesh->setModelEntity(faces[i], modelFace);
283+
// TODO: do we need to change direction of faces?
284+
}
285+
// reclassify edges
286+
Entity* f0es[3];
287+
mesh->getDownward(faces[0], 1, f0es);
288+
Entity* fes[3];
289+
mesh->getDownward(faces[1], 1, fes);
290+
for (int i = 0; i < 3; i++) {
291+
Entity* edge = fes[i];
292+
int index = findIn(f0es, 3, edge);
293+
if (index > -1) continue;
294+
mesh->setModelEntity(edge, modelFace);
295+
}
296+
mesh->getDownward(faces[2], 1, fes);
297+
for (int i = 0; i < 3; i++) {
298+
Entity* edge = fes[i];
299+
int index = findIn(f0es, 3, edge);
300+
if (index > -1) continue;
301+
mesh->setModelEntity(edge, modelFace);
302+
}
303+
Entity* vert = getTetVertOppositeFace(mesh, region, faces[0]);
304+
mesh->setModelEntity(vert, modelFace);
305+
processNewBdryVert(mesh, vert);
306+
307+
mesh->destroy(region);
308+
mesh->destroy(faces[0]);
309+
break;
310+
}
311+
case 2:
312+
{
313+
Entity* edgeToRemove;
314+
Entity* f0es[3];
315+
Entity* f1es[3];
316+
mesh->getDownward(faces[0], 1, f0es);
317+
mesh->getDownward(faces[1], 1, f1es);
318+
for (int i = 0; i < 3; i++) {
319+
edgeToRemove = f0es[i];
320+
int index = findIn(f1es, 3, edgeToRemove);
321+
if (index > -1) break;
322+
}
323+
324+
// reclassify faces
325+
for (int i = 2; i < 4; i++) {
326+
mesh->setModelEntity(faces[i], modelFace);
327+
// TODO: do we need to change direction of faces?
328+
}
329+
330+
mesh->setModelEntity(reclassifyEdge, modelFace);
331+
332+
mesh->destroy(region);
333+
mesh->destroy(faces[0]);
334+
mesh->destroy(faces[1]);
335+
mesh->destroy(edgeToRemove);
336+
break;
337+
}
338+
case 3:
339+
{
340+
mesh->setModelEntity(faces[3], modelFace);
341+
// TODO: do we need to change direction of faces?
342+
343+
Entity* vert = getTetVertOppositeFace(mesh, region, faces[3]);
344+
int n = mesh->countUpward(vert);
345+
PCU_ALWAYS_ASSERT(n == 3);
346+
Entity* edges[3];
347+
for (int i = 0; i < n; i++) {
348+
edges[i] = mesh->getUpward(vert, i);
349+
}
350+
351+
352+
mesh->destroy(region);
353+
mesh->destroy(faces[0]);
354+
mesh->destroy(faces[1]);
355+
mesh->destroy(faces[2]);
356+
mesh->destroy(edges[0]);
357+
mesh->destroy(edges[1]);
358+
mesh->destroy(edges[2]);
359+
mesh->destroy(vert);
360+
break;
361+
}
362+
}
363+
}
364+
365+
void RegionCollapse::unmark()
366+
{
367+
clearFlag(adapt,region,COLLAPSE);
368+
}
369+
370+
371+
bool setupRegionCollapse(RegionCollapse& rcollapse, Entity* region)
372+
{
373+
Adapt* adapter = rcollapse.adapt;
374+
PCU_ALWAYS_ASSERT(adapter->mesh->getType(region) == apf::Mesh::TET);
375+
if ( ! rcollapse.setRegion(region))
376+
return false;
377+
if ( ! rcollapse.checkGeom())
378+
return false;
379+
if ( ! rcollapse.checkTopo())
380+
return false;
381+
382+
return true;
383+
}
384+
385+
}

0 commit comments

Comments
 (0)