source: ThirdParty/mpqc_open/src/lib/math/isosurf/surfse.cc@ 37fe17

Adding_MD_integration_tests Adding_StructOpt_integration_tests AutomationFragmentation_failures Candidate_v1.6.1 Candidate_v1.7.0 ChemicalSpaceEvaluator Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Exclude_Hydrogens_annealWithBondGraph ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_contraction-expansion Gui_displays_atomic_force_velocity JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool PythonUI_with_named_parameters StoppableMakroAction TremoloParser_IncreasedPrecision stable
Last change on this file since 37fe17 was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 12.8 KB
RevLine 
[0b990d]1//
2// surfse.cc
3//
4// Copyright (C) 1996 Limit Point Systems, Inc.
5//
6// Author: Curtis Janssen <cljanss@limitpt.com>
7// Maintainer: LPS
8//
9// This file is part of the SC Toolkit.
10//
11// The SC Toolkit is free software; you can redistribute it and/or modify
12// it under the terms of the GNU Library General Public License as published by
13// the Free Software Foundation; either version 2, or (at your option)
14// any later version.
15//
16// The SC Toolkit is distributed in the hope that it will be useful,
17// but WITHOUT ANY WARRANTY; without even the implied warranty of
18// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19// GNU Library General Public License for more details.
20//
21// You should have received a copy of the GNU Library General Public License
22// along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24//
25// The U.S. Government is granted a limited license as per AL 91-7.
26//
27
28#include <algorithm>
29
30#include <util/misc/formio.h>
31
32#include <math/scmat/matrix.h>
33#include <math/isosurf/surf.h>
34
35using namespace std;
36using namespace sc;
37
38void
39TriangulatedSurface::remove_short_edges(double length_cutoff,
40 const Ref<Volume> &vol, double isoval)
41{
42 int j,k;
43 std::set<Ref<Triangle> >::iterator it,jt,kt;
44 std::set<Ref<Edge> >::iterator ie,je;
45
46 int surface_was_completed = _completed_surface;
47
48 if (!_completed_surface) {
49 complete_ref_arrays();
50 }
51 else {
52 clear_int_arrays();
53 }
54
55 _have_values = 0;
56 _values.clear();
57
58 if (_verbose) {
59 ExEnv::outn() << "TriangulatedSurface::remove_short_edges:" << endl
60 << "initial: ";
61 topology_info();
62 }
63
64 int deleted_edges_length;
65 do {
66 std::set<Ref<Triangle> > deleted_triangles;
67 std::set<Ref<Edge> > deleted_edges;
68 std::set<Ref<Vertex> > deleted_vertices;
69
70 std::set<Ref<Triangle> > new_triangles;
71 std::set<Ref<Edge> > new_edges;
72 std::set<Ref<Vertex> > new_vertices;
73
74 // a vertex to set-of-connected-triangles map
75 std::map<Ref<Vertex>,std::set<Ref<Triangle> > > connected_triangle_map;
76 for (it = _triangles.begin(); it != _triangles.end(); it++) {
77 Ref<Triangle> tri = *it;
78 for (j = 0; j<3; j++) {
79 Ref<Vertex> v = tri->vertex(j);
80 connected_triangle_map[v].insert(tri);
81 }
82 }
83
84 // a vertex to set-of-connected-edges map
85 std::map<Ref<Vertex>,std::set<Ref<Edge> > > connected_edge_map;
86 for (ie = _edges.begin(); ie != _edges.end(); ie++) {
87 Ref<Edge> e = *ie;
88 for (j = 0; j<2; j++) {
89 Ref<Vertex> v = e->vertex(j);
90 connected_edge_map[v].insert(e);
91 }
92 }
93
94 for (ie = _edges.begin(); ie != _edges.end(); ie++) {
95 Ref<Edge> edge = *ie;
96 double length = edge->straight_length();
97 if (length < length_cutoff) {
98 std::set<Ref<Triangle> > connected_triangles;
99 Ref<Vertex> v0 = edge->vertex(0);
100 Ref<Vertex> v1 = edge->vertex(1);
101 connected_triangles.insert(connected_triangle_map[v0].begin(),
102 connected_triangle_map[v0].end());
103 connected_triangles.insert(connected_triangle_map[v1].begin(),
104 connected_triangle_map[v1].end());
105 int skip = 0;
106 for (jt = connected_triangles.begin();
107 jt != connected_triangles.end();
108 jt++) {
109 Ref<Triangle> tri = *jt;
110 if (tri->edge(0)->straight_length() < length
111 ||tri->edge(1)->straight_length() < length
112 ||tri->edge(2)->straight_length() < length
113 ||deleted_triangles.find(tri)!=deleted_triangles.end()) {
114 skip = 1;
115 break;
116 }
117 }
118 if (skip) continue;
119 deleted_triangles.insert(connected_triangles.begin(),
120 connected_triangles.end());
121 v0 = edge->vertex(0);
122 v1 = edge->vertex(1);
123 deleted_vertices.insert(v0);
124 deleted_vertices.insert(v1);
125 // find all of the edges connected to the short edge
126 // (including the short edge)
127 std::set<Ref<Edge> > connected_edges;
128 v0 = edge->vertex(0);
129 v1 = edge->vertex(1);
130 connected_edges.insert(connected_edge_map[v0].begin(),
131 connected_edge_map[v0].end());
132 connected_edges.insert(connected_edge_map[v1].begin(),
133 connected_edge_map[v1].end());
134 deleted_edges.insert(connected_edges.begin(),
135 connected_edges.end());
136 // find the edges forming the perimeter of the deleted triangles
137 // (these are used to form the new triangles)
138 std::set<Ref<Edge> > perimeter_edges;
139 int embedded_triangle = 0;
140 for (jt = connected_triangles.begin();
141 jt != connected_triangles.end();
142 jt++) {
143 Ref<Triangle> tri = *jt;
144 for (j=0; j<3; j++) {
145 Ref<Edge> e = tri->edge(j);
146 if (connected_edges.find(e) == connected_edges.end()) {
147 // check to see if another triangle has claimed
148 // that this edge is a perimeter edge. if so
149 // then this isn't a perimeter edge after all
150 // and it must be deleted. this also implies
151 // that there is at least one triangle that
152 // has no perimeter edge. these triangles and
153 // their nonperimeter vertices must be
154 // deleted.
155 Ref<Edge> e = tri->edge(j);
156 if (perimeter_edges.find(e)
157 != perimeter_edges.end()) {
158 perimeter_edges.erase(e);
159 deleted_edges.insert(e);
160 embedded_triangle = 1;
161 }
162 else {
163 perimeter_edges.insert(e);
164 }
165 }
166 }
167 }
168 // if a triangle is embedded make sure its vertices are
169 // all deleted. (the triangle itself should already be
170 // connected and thus deleted).
171 if (embedded_triangle) {
172 // make a list of vertices on the perimeter (so i
173 // don't delete them
174 std::set<Ref<Vertex> > perimeter_vertices;
175 for (je = perimeter_edges.begin();
176 je != perimeter_edges.end();
177 je++) {
178 Ref<Edge> e = *je;
179 for (j=0; j<2; j++) {
180 Ref<Vertex> v = e->vertex(j);
181 perimeter_vertices.insert(v);
182 }
183 }
184 // find the embedded_triangle
185 for (jt = connected_triangles.begin();
186 jt != connected_triangles.end();
187 jt++) {
188 Ref<Triangle> tri = *jt;
189 // see if this triangle is embedded
190 for (j=0; j<3; j++) {
191 Ref<Edge> e = tri->edge(j);
192 if (perimeter_edges.find(e) != perimeter_edges.end())
193 break;
194 }
195 // if embedded then delete the triangle's vertices
196 if (j==3) {
197 for (j=0; j<3; j++) {
198 Ref<Vertex> v = tri->vertex(j);
199 if (perimeter_vertices.find(v) == perimeter_vertices.end())
200 deleted_vertices.insert(v);
201 }
202 }
203 }
204 }
205 // find a new point that replaces the deleted edge
206 // (for now use one of the original, since it must lie on the
207 // analytic surface)
208 Ref<Vertex> replacement_vertex = edge->vertex(0);
209 // however, if we have a volume, find a new vertex on
210 // the analytic surface near the center of the edge
211 if (vol.nonnull()) {
212 SCVector3 point, norm;
213 int hn = edge->interpolate(0.5,point,norm,vol,isoval);
214 replacement_vertex = new Vertex(point);
215 if (hn) replacement_vertex->set_normal(norm);
216 }
217 new_vertices.insert(replacement_vertex);
218 // for each vertex on the perimeter form a new edge to the
219 // replacement vertex
220 std::map<Ref<Vertex>,Ref<Edge> > new_edge_map;
221 for (je = perimeter_edges.begin(); je!=perimeter_edges.end();
222 je++) {
223 Ref<Edge> e = *je;
224 for (k = 0; k<2; k++) {
225 Ref<Vertex> v = e->vertex(k);
226 if (new_edge_map.find(v) == new_edge_map.end()) {
227 Ref<Edge> new_e = newEdge(
228 replacement_vertex,
229 v
230 );
231 new_edge_map[v] = new_e;
232 new_edges.insert(new_e);
233 }
234 }
235 }
236 // for each edge on the perimeter form a new triangle with the
237 // replacement vertex
238 for (je = perimeter_edges.begin(); je != perimeter_edges.end();
239 je++) {
240 Ref<Edge> e1 = *je;
241 Ref<Vertex> v0 = e1->vertex(0);
242 Ref<Vertex> v1 = e1->vertex(1);
243 Ref<Edge> e2 = new_edge_map[v0];
244 Ref<Edge> e3 = new_edge_map[v1];
245 if (e1.null() || e2.null() || e3.null()) {
246 ExEnv::errn() << "TriangulatedSurface::remove_short_edges: "
247 << "building new triangle but edges are null:"
248 << endl;
249 if (e1.null()) ExEnv::errn() << " e1" << endl;
250 if (e2.null()) ExEnv::errn() << " e2" << endl;
251 if (e3.null()) ExEnv::errn() << " e3" << endl;
252 abort();
253 }
254 // Compute the correct orientation of e1 within the new
255 // triangle, by finding the orientation within the old
256 // triangle.
257 int orientation = 0;
258 for (kt = connected_triangles.begin();
259 kt != connected_triangles.end();
260 kt++) {
261 if ((*kt)->contains(e1)) {
262 orientation
263 = (*kt)->orientation(e1);
264 break;
265 }
266 }
267 Ref<Triangle> newtri(newTriangle(e1,e2,e3,orientation));
268 new_triangles.insert(newtri);
269 }
270 //ExEnv::outn() << "WARNING: only one short edge removed" << endl;
271 //break;
272 }
273 }
274
275 erase_elements_by_value(_triangles,
276 deleted_triangles.begin(),
277 deleted_triangles.end());
278
279 erase_elements_by_value(_edges,
280 deleted_edges.begin(),
281 deleted_edges.end());
282
283 erase_elements_by_value(_vertices,
284 deleted_vertices.begin(),
285 deleted_vertices.end());
286
287 _triangles.insert(new_triangles.begin(), new_triangles.end());
288 _edges.insert(new_edges.begin(), new_edges.end());
289 _vertices.insert(new_vertices.begin(), new_vertices.end());
290
291 if (_verbose) {
292 topology_info();
293 }
294
295 deleted_edges_length = deleted_edges.size();
296 //ExEnv::outn() << "WARNING: one pass short edge removal" << endl;
297 //deleted_edges_length = 0; // do one pass
298 } while(deleted_edges_length != 0);
299
300 // fix the index maps
301 recompute_index_maps();
302
303 // restore the int arrays if they were there to begin with
304 if (surface_was_completed) {
305 complete_int_arrays();
306 _completed_surface = 1;
307 }
308
309 if (_verbose) {
310 ExEnv::outn() << "final: ";
311 topology_info();
312 }
313}
314
315/////////////////////////////////////////////////////////////////////////////
316
317// Local Variables:
318// mode: c++
319// c-file-style: "CLJ"
320// End:
Note: See TracBrowser for help on using the repository browser.