ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/applications/nanoRodBuilder/GeometryBuilder.cpp
Revision: 2218
Committed: Tue May 3 17:55:28 2005 UTC (19 years, 2 months ago) by chuckv
File size: 11889 byte(s)
Log Message:
More changes to nanoRodBuilder....

File Contents

# User Rev Content
1 chuckv 2164 /*
2     * GeometryBuilder.cpp
3     * nanorodBuilder
4     *
5     * Created by Charles Vardeman II on 4/4/05.
6     * Copyright 2005 University of Notre Dame. All rights reserved.
7     *
8     */
9    
10     /*
11     * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
12     *
13     * The University of Notre Dame grants you ("Licensee") a
14     * non-exclusive, royalty free, license to use, modify and
15     * redistribute this software in source and binary code form, provided
16     * that the following conditions are met:
17     *
18     * 1. Acknowledgement of the program authors must be made in any
19     * publication of scientific results based in part on use of the
20     * program. An acceptable form of acknowledgement is citation of
21     * the article in which the program was described (Matthew
22 gezelter 2204 * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
23     * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
24 chuckv 2164 * Parallel Simulation Engine for Molecular Dynamics,"
25 gezelter 2204 * J. Comput. Chem. 26, pp. 252-271 (2005))
26 chuckv 2164 *
27     * 2. Redistributions of source code must retain the above copyright
28     * notice, this list of conditions and the following disclaimer.
29     *
30     * 3. Redistributions in binary form must reproduce the above copyright
31     * notice, this list of conditions and the following disclaimer in the
32     * documentation and/or other materials provided with the
33     * distribution.
34     *
35     * This software is provided "AS IS," without a warranty of any
36     * kind. All express or implied conditions, representations and
37     * warranties, including any implied warranty of merchantability,
38     * fitness for a particular purpose or non-infringement, are hereby
39     * excluded. The University of Notre Dame and its licensors shall not
40     * be liable for any damages suffered by licensee as a result of
41     * using, modifying or distributing the software or its
42     * derivatives. In no event will the University of Notre Dame or its
43     * licensors be liable for any lost revenue, profit or data, or for
44     * direct, indirect, special, consequential, incidental or punitive
45     * damages, however caused and regardless of the theory of liability,
46     * arising out of the use of or inability to use software, even if the
47     * University of Notre Dame has been advised of the possibility of
48     * such damages.
49     */
50    
51     #include "GeometryBuilder.hpp"
52    
53     #include <CGAL/Simple_cartesian.h>
54     #include <CGAL/Polyhedron_incremental_builder_3.h>
55     #include <CGAL/Polyhedron_3.h>
56     #include <CGAL/Homogeneous.h>
57     #include <CGAL/Polyhedron_traits_with_normals_3.h>
58     #include <CGAL/Polyhedron_3.h>
59     #include <CGAL/Aff_transformation_3.h>
60     #include <CGAL/aff_transformation_tags.h>
61     #include <iostream>
62     #include <algorithm>
63 chuckv 2215 #include <fstream>
64     #include <math.h>
65 chuckv 2164
66 chuckv 2215 using namespace std;
67 chuckv 2164 using namespace oopse;
68    
69     //typedef CGAL::Homogeneous<int> Kernel;
70     typedef CGAL::Simple_cartesian<double> Kernel;
71     //typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
72    
73     typedef Kernel::Point_3 Point_3;
74     typedef Kernel::Vector_3 Vector_3;
75     typedef CGAL::Polyhedron_traits_with_normals_3<Kernel> Traits;
76     //typedef CGAL::Polyhedron_3<Kernel> Polyhedron;
77     typedef CGAL::Polyhedron_3<Traits> Polyhedron;
78     typedef Polyhedron::HalfedgeDS HalfedgeDS;
79     typedef Polyhedron::Facet_iterator Facet_iterator;
80     typedef Polyhedron::Halfedge_around_facet_circulator Halfedge_facet_circulator;
81     typedef Polyhedron::Facet_iterator Facet_iterator;
82     typedef Polyhedron::Plane_iterator Plane_iterator;
83     typedef Polyhedron::Vertex_handle Vertex_handle;
84    
85     Polyhedron nanoRodPolyhedron;
86    
87    
88    
89    
90    
91    
92    
93     //typedef CGAL::Scaling Scaling;
94     //typedef Aff_transformation_3<Kernel> A;( const Scaling,
95     // Kernel::RT s=RT(20),
96     // Kernel::RT hw = RT(1));
97    
98     // A modifier creating a triangle with the incremental builder.
99     template <class HDS>
100 chuckv 2218 class buildSingleCrystal : public CGAL::Modifier_base<HDS> {
101 gezelter 2204 public:
102 chuckv 2164 Vertex_handle end1;
103     Vertex_handle neight1;
104 gezelter 2204 Vertex_handle end2;
105 chuckv 2164 Vertex_handle neight2;
106 chuckv 2215 Vertex_handle neight3;
107 chuckv 2164
108 chuckv 2218 buildSingleCrystal() {}
109 chuckv 2164 void operator()( HDS& hds) {
110     // Postcondition: `hds' is a valid polyhedral surface.
111     CGAL::Polyhedron_incremental_builder_3<HDS> B( hds, true);
112     B.begin_surface( 12, 15, 6);
113     typedef typename HDS::Vertex Vertex;
114     typedef typename Vertex::Point Point;
115    
116    
117    
118    
119    
120     B.add_vertex( Point(-0.7887222926324, 0.4874571845315, -0.2562714077342));
121     B.add_vertex( Point(-0.4874571845316, 0.4874571845315, 0.6709272557930));
122     B.add_vertex( Point(-0.7887222926324, -0.4874571845315, -0.2562714077342)); //End vertex
123     end1 = B.add_vertex( Point( 0.0000000000000, 1.0000000000000, 0.0000000000000));
124 chuckv 2215 neight3 = B.add_vertex( Point(-0.4874571845315, -0.4874571845316, 0.6709272557930));
125 chuckv 2164 neight1 = B.add_vertex( Point(-0.0000000000000, 0.4874571845316, -0.8293116961175));
126     B.add_vertex( Point( 0.0000000000000, -0.4874571845316, -0.8293116961175));
127     B.add_vertex( Point( 0.4874571845315, 0.4874571845316, 0.6709272557930));
128     end2 = B.add_vertex( Point(-0.0000000000000, -1.0000000000000, 0.0000000000000)); //End Vertex
129     B.add_vertex( Point( 0.7887222926324, 0.4874571845315, -0.2562714077342));
130     neight2 = B.add_vertex( Point( 0.4874571845316, -0.4874571845315, 0.6709272557930));
131     B.add_vertex( Point( 0.7887222926324, -0.4874571845315, -0.2562714077342));
132    
133     B.begin_facet();
134     B.add_vertex_to_facet( 7);
135     B.add_vertex_to_facet( 9);
136     B.add_vertex_to_facet( 11);
137     B.add_vertex_to_facet( 10);
138     B.end_facet();
139    
140     B.begin_facet();
141     B.add_vertex_to_facet( 8);
142     B.add_vertex_to_facet( 10);
143     B.add_vertex_to_facet( 11);
144     B.end_facet();
145    
146     B.begin_facet();
147     B.add_vertex_to_facet( 3);
148     B.add_vertex_to_facet( 9);
149     B.add_vertex_to_facet( 7);
150     B.end_facet();
151    
152     B.begin_facet();
153     B.add_vertex_to_facet( 9);
154     B.add_vertex_to_facet( 5);
155     B.add_vertex_to_facet( 6);
156     B.add_vertex_to_facet( 11);
157     B.end_facet();
158    
159     B.begin_facet();
160     B.add_vertex_to_facet( 8);
161     B.add_vertex_to_facet( 11);
162     B.add_vertex_to_facet( 6);
163     B.end_facet();
164    
165     B.begin_facet();
166     B.add_vertex_to_facet( 3);
167     B.add_vertex_to_facet( 5);
168     B.add_vertex_to_facet( 9);
169     B.end_facet();
170    
171     B.begin_facet();
172     B.add_vertex_to_facet( 5);
173     B.add_vertex_to_facet( 0);
174     B.add_vertex_to_facet( 2);
175     B.add_vertex_to_facet( 6);
176     B.end_facet();
177    
178     B.begin_facet();
179     B.add_vertex_to_facet( 8);
180     B.add_vertex_to_facet( 6);
181     B.add_vertex_to_facet( 2);
182     B.end_facet();
183    
184     B.begin_facet();
185     B.add_vertex_to_facet( 3);
186     B.add_vertex_to_facet( 0);
187     B.add_vertex_to_facet( 5);
188     B.end_facet();
189    
190     B.begin_facet();
191     B.add_vertex_to_facet( 0);
192     B.add_vertex_to_facet( 1);
193     B.add_vertex_to_facet( 4);
194     B.add_vertex_to_facet( 2);
195     B.end_facet();
196    
197     B.begin_facet();
198     B.add_vertex_to_facet( 8);
199     B.add_vertex_to_facet( 2);
200     B.add_vertex_to_facet( 4);
201     B.end_facet();
202    
203     B.begin_facet();
204     B.add_vertex_to_facet( 3);
205     B.add_vertex_to_facet( 1);
206     B.add_vertex_to_facet( 0);
207     B.end_facet();
208    
209     B.begin_facet();
210     B.add_vertex_to_facet( 1);
211     B.add_vertex_to_facet( 7);
212     B.add_vertex_to_facet( 10);
213     B.add_vertex_to_facet( 4);
214     B.end_facet();
215    
216     B.begin_facet();
217     B.add_vertex_to_facet( 8);
218     B.add_vertex_to_facet( 4);
219     B.add_vertex_to_facet( 10);
220     B.end_facet();
221    
222     B.begin_facet();
223     B.add_vertex_to_facet( 3);
224     B.add_vertex_to_facet( 7);
225     B.add_vertex_to_facet( 1);
226     B.end_facet();
227    
228    
229     B.end_surface();
230     }
231     };
232    
233    
234    
235     struct Normal_vector {
236     template <class Facet>
237     typename Facet::Plane_3 operator()( Facet& f) {
238     typename Facet::Halfedge_handle h = f.halfedge();
239     // Facet::Plane_3 is the normal vector type. We assume the
240     // CGAL Kernel here and use its global functions.
241     return CGAL::cross_product(
242     h->next()->vertex()->point() - h->vertex()->point(),
243     h->next()->next()->vertex()->point() - h->next()->vertex()->point());
244     }
245     };
246    
247    
248     bool GeometryBuilder::isInsidePolyhedron(double x, double y, double z) {
249    
250 gezelter 2204 Point_3 point(x,y,z);
251 chuckv 2164 Plane_iterator i;
252     Facet_iterator j;
253     for ( i =nanoRodPolyhedron.planes_begin(), j = nanoRodPolyhedron.facets_begin(); i != nanoRodPolyhedron.planes_end() && j !=nanoRodPolyhedron.facets_end(); ++i, ++j) {
254     Halfedge_facet_circulator k = j->facet_begin();
255    
256     Vector_3 newVector = point - k->vertex()->point();
257     Vector_3 normal = *i;
258     double dot_product = newVector.x() * normal.x() + newVector.y() * normal.y() + newVector.z() * normal.z();
259    
260     if (dot_product < 0) {
261     return false;
262     }
263     }
264    
265     return true;
266     }
267    
268    
269     GeometryBuilder::GeometryBuilder(double length,double width) {
270 gezelter 2204 // Create the geometry for nanorod
271 chuckv 2218 buildSingleCrystal<HalfedgeDS> nanorod;
272 chuckv 2164
273     nanoRodPolyhedron.delegate( nanorod);
274    
275 gezelter 2204 double y1 = nanorod.end1->point().y() - nanorod.neight1->point().y();
276     double y2 = nanorod.end2->point().y() - nanorod.neight2->point().y();
277 chuckv 2215
278     double endDist = sqrt(pow(nanorod.neight2->point().x() - nanorod.neight3->point().x(),2)+
279     pow(nanorod.neight2->point().y() - nanorod.neight3->point().y(),2)+
280     pow(nanorod.neight2->point().z() - nanorod.neight3->point().z(),2));
281    
282     double endRatio1 = y1/endDist;
283     double endRatio2 = y2/endDist;
284    
285     std::cout << "End dist is " << endDist <<" ratio " << endRatio1 << std::endl;
286    
287 gezelter 2204 CGAL::Aff_transformation_3<Kernel> aff_tranformation( width,
288     0.0,
289     0.0,
290     0.0,
291     0.0,
292     length,
293     0.0,
294     0.0,
295     0.0,
296     0.0,
297     width,
298     0.0);
299     std::transform( nanoRodPolyhedron.points_begin(), nanoRodPolyhedron.points_end(), nanoRodPolyhedron.points_begin(), aff_tranformation);
300 chuckv 2164
301 chuckv 2215
302     double endDist2 = sqrt(pow(nanorod.neight2->point().x() - nanorod.neight3->point().x(),2)+
303     pow(nanorod.neight2->point().y() - nanorod.neight3->point().y(),2)+
304     pow(nanorod.neight2->point().z() - nanorod.neight3->point().z(),2));
305    
306     Point_3 point1(nanorod.end1->point().x(), endDist2*endRatio1 + nanorod.neight1->point().y(), nanorod.end1->point().z());
307     Point_3 point2(nanorod.end2->point().x(), endDist2*endRatio2 + nanorod.neight2->point().y(), nanorod.end2->point().z());
308 gezelter 2204 nanorod.end1->point() = point1;
309     nanorod.end2->point() = point2;
310 chuckv 2164
311 gezelter 2204 // Construct normal vector for each face.
312     std::transform( nanoRodPolyhedron.facets_begin(), nanoRodPolyhedron.facets_end(), nanoRodPolyhedron.planes_begin(),
313     Normal_vector());
314 chuckv 2215 }
315     void GeometryBuilder::dumpGeometry(const std::string& geomFileName){
316    
317     std::ofstream newGeomFile;
318    
319     //create new .md file based on old .md file
320     newGeomFile.open(geomFileName.c_str());
321    
322     // Write polyhedron in Object File Format (OFF).
323     CGAL::set_ascii_mode( std::cout);
324     newGeomFile << "OFF" << std::endl << nanoRodPolyhedron.size_of_vertices() << ' '
325     << nanoRodPolyhedron.size_of_facets() << " 0" << std::endl;
326     std::copy( nanoRodPolyhedron.points_begin(), nanoRodPolyhedron.points_end(),
327     std::ostream_iterator<Point_3>( newGeomFile, "\n"));
328     for ( Facet_iterator i = nanoRodPolyhedron.facets_begin(); i != nanoRodPolyhedron.facets_end(); ++i) {
329     Halfedge_facet_circulator j = i->facet_begin();
330     // Facets in polyhedral surfaces are at least triangles.
331     CGAL_assertion( CGAL::circulator_size(j) >= 3);
332     newGeomFile << CGAL::circulator_size(j) << ' ';
333     do {
334     newGeomFile << ' ' << std::distance(nanoRodPolyhedron.vertices_begin(), j->vertex());
335     } while ( ++j != i->facet_begin());
336     newGeomFile << std::endl;
337     }
338    
339     newGeomFile.close();
340    
341    
342 chuckv 2164
343    
344 chuckv 2215 }
345 chuckv 2164
346 chuckv 2215