applications/nanoRodBuilder/GeometryBuilder.cpp

/*
 *  GeometryBuilder.cpp
 *  nanorodBuilder
 *
 *  Created by Charles Vardeman II on 4/4/05.
 *  Copyright 2005 University of Notre Dame. All rights reserved.
 *
 */

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */

#include "GeometryBuilder.hpp"

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Polyhedron_incremental_builder_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Homogeneous.h>
#include <CGAL/Polyhedron_traits_with_normals_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Aff_transformation_3.h>
#include <CGAL/aff_transformation_tags.h>
#include <iostream>
#include <algorithm>
#include <fstream>
#include <math.h>

using namespace std;
using namespace oopse;

//typedef CGAL::Homogeneous<int>              Kernel;
typedef CGAL::Simple_cartesian<double>     Kernel;
//typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;

typedef Kernel::Point_3                    Point_3;
typedef Kernel::Vector_3                               Vector_3;
typedef CGAL::Polyhedron_traits_with_normals_3<Kernel> Traits;
//typedef CGAL::Polyhedron_3<Kernel>                     Polyhedron;
typedef CGAL::Polyhedron_3<Traits>                     Polyhedron;
typedef Polyhedron::HalfedgeDS             HalfedgeDS;
typedef Polyhedron::Facet_iterator                   Facet_iterator;
typedef Polyhedron::Halfedge_around_facet_circulator Halfedge_facet_circulator;
typedef Polyhedron::Facet_iterator Facet_iterator;
typedef Polyhedron::Plane_iterator Plane_iterator;
typedef Polyhedron::Vertex_handle Vertex_handle;

Polyhedron nanoRodPolyhedron;


//typedef CGAL::Scaling Scaling;
//typedef Aff_transformation_3<Kernel> A;( const Scaling,
//                                                                                                                              Kernel::RT s=RT(20),
//                                                                                                                              Kernel::RT hw = RT(1));


// A modifier creating a triangle with the incremental builder.
template <class HDS>
class buildSingleCrystal : public CGAL::Modifier_base<HDS> {
 public:
  Vertex_handle end1;
  Vertex_handle neight1;
  Vertex_handle end2;
  Vertex_handle neight2;
  Vertex_handle neight3;
  
  buildSingleCrystal() {}
  void operator()( HDS& hds) {
    // Postcondition: `hds' is a valid polyhedral surface.
    CGAL::Polyhedron_incremental_builder_3<HDS> B( hds, true);
    B.begin_surface( 12, 15, 6);
    typedef typename HDS::Vertex   Vertex;
    typedef typename Vertex::Point Point;
    
    
    B.add_vertex( Point(-0.7887222926324,    0.4874571845315, -0.2562714077342));
    B.add_vertex( Point(-0.4874571845316,    0.4874571845315,  0.6709272557930));
    B.add_vertex( Point(-0.7887222926324,   -0.4874571845315, -0.2562714077342)); //End vertex
    end1 = B.add_vertex( Point( 0.0000000000000,    1.0000000000000,  0.0000000000000));
    neight3 = B.add_vertex( Point(-0.4874571845315,    -0.4874571845316,  0.6709272557930));
    neight1 = B.add_vertex( Point(-0.0000000000000,    0.4874571845316, -0.8293116961175));
    B.add_vertex( Point( 0.0000000000000,    -0.4874571845316, -0.8293116961175));
    B.add_vertex( Point( 0.4874571845315,    0.4874571845316,  0.6709272557930));
    end2 = B.add_vertex( Point(-0.0000000000000,    -1.0000000000000,  0.0000000000000)); //End Vertex
    B.add_vertex( Point( 0.7887222926324,    0.4874571845315, -0.2562714077342));
    neight2 = B.add_vertex( Point( 0.4874571845316,    -0.4874571845315,  0.6709272557930));
    B.add_vertex( Point( 0.7887222926324,    -0.4874571845315, -0.2562714077342));
    
    B.begin_facet();
    B.add_vertex_to_facet( 7);
    B.add_vertex_to_facet( 9);
    B.add_vertex_to_facet( 11);
    B.add_vertex_to_facet( 10);
    B.end_facet();
                                
    B.begin_facet();
    B.add_vertex_to_facet( 8);
    B.add_vertex_to_facet( 10);
    B.add_vertex_to_facet( 11);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 3);
    B.add_vertex_to_facet( 9);
    B.add_vertex_to_facet( 7);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 9);
    B.add_vertex_to_facet( 5);
    B.add_vertex_to_facet( 6);
    B.add_vertex_to_facet( 11);
    B.end_facet();
                                
    B.begin_facet();
    B.add_vertex_to_facet( 8);
    B.add_vertex_to_facet( 11);
    B.add_vertex_to_facet( 6);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 3);
    B.add_vertex_to_facet( 5);
    B.add_vertex_to_facet( 9);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 5);
    B.add_vertex_to_facet( 0);
    B.add_vertex_to_facet( 2);
    B.add_vertex_to_facet( 6);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 8);
    B.add_vertex_to_facet( 6);
    B.add_vertex_to_facet( 2);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 3);
    B.add_vertex_to_facet( 0);
    B.add_vertex_to_facet( 5);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 0);
    B.add_vertex_to_facet( 1);
    B.add_vertex_to_facet( 4);
    B.add_vertex_to_facet( 2);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 8);
    B.add_vertex_to_facet( 2);
    B.add_vertex_to_facet( 4);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 3);
    B.add_vertex_to_facet( 1);
    B.add_vertex_to_facet( 0);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 1);
    B.add_vertex_to_facet( 7);
    B.add_vertex_to_facet( 10);
    B.add_vertex_to_facet( 4);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 8);
    B.add_vertex_to_facet( 4);
    B.add_vertex_to_facet( 10);
    B.end_facet();
    
    B.begin_facet();
    B.add_vertex_to_facet( 3);
    B.add_vertex_to_facet( 7);
    B.add_vertex_to_facet( 1);
    B.end_facet();
    
                                
    B.end_surface();
  }
};


struct Normal_vector {
  template <class Facet>
  typename Facet::Plane_3 operator()( Facet& f) {
    typename Facet::Halfedge_handle h = f.halfedge();
    // Facet::Plane_3 is the normal vector type. We assume the
    // CGAL Kernel here and use its global functions.
    return CGAL::cross_product( 
                               h->next()->vertex()->point() - h->vertex()->point(),
                               h->next()->next()->vertex()->point() - h->next()->vertex()->point());
  }
};


bool GeometryBuilder::isInsidePolyhedron(double x, double y, double z) {
        
  Point_3 point(x,y,z);
  Plane_iterator i;
  Facet_iterator j;
  for ( i =nanoRodPolyhedron.planes_begin(), j = nanoRodPolyhedron.facets_begin(); i != nanoRodPolyhedron.planes_end() && j !=nanoRodPolyhedron.facets_end(); ++i, ++j) {
    Halfedge_facet_circulator k = j->facet_begin();
    
    Vector_3 newVector = point - k->vertex()->point();
    Vector_3 normal = *i;               
    double dot_product = newVector.x() * normal.x() + newVector.y() * normal.y() + newVector.z() * normal.z();
    
    if (dot_product < 0) {
      return false;     
    }
  }
  
  return true;
}


GeometryBuilder::GeometryBuilder(double length,double width) {
  // Create the geometry for nanorod
 buildSingleCrystal<HalfedgeDS> nanorod;
  
  nanoRodPolyhedron.delegate( nanorod);
   
  double y1 = nanorod.end1->point().y() - nanorod.neight1->point().y();
  double y2 = nanorod.end2->point().y() - nanorod.neight2->point().y();
  
  double endDist = sqrt(pow(nanorod.neight2->point().x() - nanorod.neight3->point().x(),2)+
                        pow(nanorod.neight2->point().y() - nanorod.neight3->point().y(),2)+
                        pow(nanorod.neight2->point().z() - nanorod.neight3->point().z(),2));
  
  double endRatio1 = y1/endDist;
  double endRatio2 = y2/endDist;
  
  std::cout << "End dist is " << endDist <<" ratio " << endRatio1 << std::endl;
  
  CGAL::Aff_transformation_3<Kernel> aff_tranformation( width,
                                                        0.0,
                                                        0.0,
                                                        0.0,
                                                        0.0,
                                                        length,
                                                        0.0,
                                                        0.0,
                                                        0.0,
                                                        0.0,
                                                        width,
                                                        0.0);   
  std::transform( nanoRodPolyhedron.points_begin(), nanoRodPolyhedron.points_end(), nanoRodPolyhedron.points_begin(), aff_tranformation);
        
  
  double endDist2 = sqrt(pow(nanorod.neight2->point().x() - nanorod.neight3->point().x(),2)+
                        pow(nanorod.neight2->point().y() - nanorod.neight3->point().y(),2)+
                        pow(nanorod.neight2->point().z() - nanorod.neight3->point().z(),2));
  
  Point_3 point1(nanorod.end1->point().x(), endDist2*endRatio1 + nanorod.neight1->point().y(), nanorod.end1->point().z());
  Point_3 point2(nanorod.end2->point().x(), endDist2*endRatio2 + nanorod.neight2->point().y(), nanorod.end2->point().z());
  nanorod.end1->point() = point1;
  nanorod.end2->point() = point2;
        
  // Construct normal vector for each face.
  std::transform( nanoRodPolyhedron.facets_begin(), nanoRodPolyhedron.facets_end(), nanoRodPolyhedron.planes_begin(),
                  Normal_vector());
}       
  void GeometryBuilder::dumpGeometry(const std::string& geomFileName){
     
     std::ofstream newGeomFile;
     
     //create new .md file based on old .md file
     newGeomFile.open(geomFileName.c_str());
     
     // Write polyhedron in Object File Format (OFF).
     CGAL::set_ascii_mode( std::cout);
     newGeomFile << "OFF" << std::endl << nanoRodPolyhedron.size_of_vertices() << ' ' 
        << nanoRodPolyhedron.size_of_facets() << " 0" << std::endl;
     std::copy( nanoRodPolyhedron.points_begin(), nanoRodPolyhedron.points_end(),
                std::ostream_iterator<Point_3>( newGeomFile, "\n"));
     for (  Facet_iterator i = nanoRodPolyhedron.facets_begin(); i != nanoRodPolyhedron.facets_end(); ++i) {
        Halfedge_facet_circulator j = i->facet_begin();
        // Facets in polyhedral surfaces are at least triangles.
        CGAL_assertion( CGAL::circulator_size(j) >= 3);
        newGeomFile << CGAL::circulator_size(j) << ' ';
        do {
           newGeomFile << ' ' << std::distance(nanoRodPolyhedron.vertices_begin(), j->vertex());
        } while ( ++j != i->facet_begin());
        newGeomFile << std::endl;
     }
     
     newGeomFile.close();
     
     
  }


Revision:	2239
Committed:	Sun May 22 21:36:21 2005 UTC (19 years, 1 month ago) by chuckv
File size:	11893 byte(s)
Log Message:	Fix to print pressure tensor.
#	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	chuckv	2239
99
100
101
102	chuckv	2164	// A modifier creating a triangle with the incremental builder.
103			template <class HDS>
104	chuckv	2218	class buildSingleCrystal : public CGAL::Modifier_base<HDS> {
105	gezelter	2204	public:
106	chuckv	2164	Vertex_handle end1;
107			Vertex_handle neight1;
108	gezelter	2204	Vertex_handle end2;
109	chuckv	2164	Vertex_handle neight2;
110	chuckv	2215	Vertex_handle neight3;
111	chuckv	2164
112	chuckv	2218	buildSingleCrystal() {}
113	chuckv	2164	void operator()( HDS& hds) {
114			// Postcondition: `hds' is a valid polyhedral surface.
115			CGAL::Polyhedron_incremental_builder_3<HDS> B( hds, true);
116			B.begin_surface( 12, 15, 6);
117			typedef typename HDS::Vertex Vertex;
118			typedef typename Vertex::Point Point;
119
120
121
122
123
124			B.add_vertex( Point(-0.7887222926324, 0.4874571845315, -0.2562714077342));
125			B.add_vertex( Point(-0.4874571845316, 0.4874571845315, 0.6709272557930));
126			B.add_vertex( Point(-0.7887222926324, -0.4874571845315, -0.2562714077342)); //End vertex
127			end1 = B.add_vertex( Point( 0.0000000000000, 1.0000000000000, 0.0000000000000));
128	chuckv	2215	neight3 = B.add_vertex( Point(-0.4874571845315, -0.4874571845316, 0.6709272557930));
129	chuckv	2164	neight1 = B.add_vertex( Point(-0.0000000000000, 0.4874571845316, -0.8293116961175));
130			B.add_vertex( Point( 0.0000000000000, -0.4874571845316, -0.8293116961175));
131			B.add_vertex( Point( 0.4874571845315, 0.4874571845316, 0.6709272557930));
132			end2 = B.add_vertex( Point(-0.0000000000000, -1.0000000000000, 0.0000000000000)); //End Vertex
133			B.add_vertex( Point( 0.7887222926324, 0.4874571845315, -0.2562714077342));
134			neight2 = B.add_vertex( Point( 0.4874571845316, -0.4874571845315, 0.6709272557930));
135			B.add_vertex( Point( 0.7887222926324, -0.4874571845315, -0.2562714077342));
136
137			B.begin_facet();
138			B.add_vertex_to_facet( 7);
139			B.add_vertex_to_facet( 9);
140			B.add_vertex_to_facet( 11);
141			B.add_vertex_to_facet( 10);
142			B.end_facet();
143
144			B.begin_facet();
145			B.add_vertex_to_facet( 8);
146			B.add_vertex_to_facet( 10);
147			B.add_vertex_to_facet( 11);
148			B.end_facet();
149
150			B.begin_facet();
151			B.add_vertex_to_facet( 3);
152			B.add_vertex_to_facet( 9);
153			B.add_vertex_to_facet( 7);
154			B.end_facet();
155
156			B.begin_facet();
157			B.add_vertex_to_facet( 9);
158			B.add_vertex_to_facet( 5);
159			B.add_vertex_to_facet( 6);
160			B.add_vertex_to_facet( 11);
161			B.end_facet();
162
163			B.begin_facet();
164			B.add_vertex_to_facet( 8);
165			B.add_vertex_to_facet( 11);
166			B.add_vertex_to_facet( 6);
167			B.end_facet();
168
169			B.begin_facet();
170			B.add_vertex_to_facet( 3);
171			B.add_vertex_to_facet( 5);
172			B.add_vertex_to_facet( 9);
173			B.end_facet();
174
175			B.begin_facet();
176			B.add_vertex_to_facet( 5);
177			B.add_vertex_to_facet( 0);
178			B.add_vertex_to_facet( 2);
179			B.add_vertex_to_facet( 6);
180			B.end_facet();
181
182			B.begin_facet();
183			B.add_vertex_to_facet( 8);
184			B.add_vertex_to_facet( 6);
185			B.add_vertex_to_facet( 2);
186			B.end_facet();
187
188			B.begin_facet();
189			B.add_vertex_to_facet( 3);
190			B.add_vertex_to_facet( 0);
191			B.add_vertex_to_facet( 5);
192			B.end_facet();
193
194			B.begin_facet();
195			B.add_vertex_to_facet( 0);
196			B.add_vertex_to_facet( 1);
197			B.add_vertex_to_facet( 4);
198			B.add_vertex_to_facet( 2);
199			B.end_facet();
200
201			B.begin_facet();
202			B.add_vertex_to_facet( 8);
203			B.add_vertex_to_facet( 2);
204			B.add_vertex_to_facet( 4);
205			B.end_facet();
206
207			B.begin_facet();
208			B.add_vertex_to_facet( 3);
209			B.add_vertex_to_facet( 1);
210			B.add_vertex_to_facet( 0);
211			B.end_facet();
212
213			B.begin_facet();
214			B.add_vertex_to_facet( 1);
215			B.add_vertex_to_facet( 7);
216			B.add_vertex_to_facet( 10);
217			B.add_vertex_to_facet( 4);
218			B.end_facet();
219
220			B.begin_facet();
221			B.add_vertex_to_facet( 8);
222			B.add_vertex_to_facet( 4);
223			B.add_vertex_to_facet( 10);
224			B.end_facet();
225
226			B.begin_facet();
227			B.add_vertex_to_facet( 3);
228			B.add_vertex_to_facet( 7);
229			B.add_vertex_to_facet( 1);
230			B.end_facet();
231
232
233			B.end_surface();
234			}
235			};
236
237
238
239			struct Normal_vector {
240			template <class Facet>
241			typename Facet::Plane_3 operator()( Facet& f) {
242			typename Facet::Halfedge_handle h = f.halfedge();
243			// Facet::Plane_3 is the normal vector type. We assume the
244			// CGAL Kernel here and use its global functions.
245			return CGAL::cross_product(
246			h->next()->vertex()->point() - h->vertex()->point(),
247			h->next()->next()->vertex()->point() - h->next()->vertex()->point());
248			}
249			};
250
251
252			bool GeometryBuilder::isInsidePolyhedron(double x, double y, double z) {
253
254	gezelter	2204	Point_3 point(x,y,z);
255	chuckv	2164	Plane_iterator i;
256			Facet_iterator j;
257			for ( i =nanoRodPolyhedron.planes_begin(), j = nanoRodPolyhedron.facets_begin(); i != nanoRodPolyhedron.planes_end() && j !=nanoRodPolyhedron.facets_end(); ++i, ++j) {
258			Halfedge_facet_circulator k = j->facet_begin();
259
260			Vector_3 newVector = point - k->vertex()->point();
261			Vector_3 normal = *i;
262			double dot_product = newVector.x() * normal.x() + newVector.y() * normal.y() + newVector.z() * normal.z();
263
264			if (dot_product < 0) {
265			return false;
266			}
267			}
268
269			return true;
270			}
271
272
273			GeometryBuilder::GeometryBuilder(double length,double width) {
274	gezelter	2204	// Create the geometry for nanorod
275	chuckv	2218	buildSingleCrystal<HalfedgeDS> nanorod;
276	chuckv	2164
277			nanoRodPolyhedron.delegate( nanorod);
278
279	gezelter	2204	double y1 = nanorod.end1->point().y() - nanorod.neight1->point().y();
280			double y2 = nanorod.end2->point().y() - nanorod.neight2->point().y();
281	chuckv	2215
282			double endDist = sqrt(pow(nanorod.neight2->point().x() - nanorod.neight3->point().x(),2)+
283			pow(nanorod.neight2->point().y() - nanorod.neight3->point().y(),2)+
284			pow(nanorod.neight2->point().z() - nanorod.neight3->point().z(),2));
285
286			double endRatio1 = y1/endDist;
287			double endRatio2 = y2/endDist;
288
289			std::cout << "End dist is " << endDist <<" ratio " << endRatio1 << std::endl;
290
291	gezelter	2204	CGAL::Aff_transformation_3<Kernel> aff_tranformation( width,
292			0.0,
293			0.0,
294			0.0,
295			0.0,
296			length,
297			0.0,
298			0.0,
299			0.0,
300			0.0,
301			width,
302			0.0);
303			std::transform( nanoRodPolyhedron.points_begin(), nanoRodPolyhedron.points_end(), nanoRodPolyhedron.points_begin(), aff_tranformation);
304	chuckv	2164
305	chuckv	2215
306			double endDist2 = sqrt(pow(nanorod.neight2->point().x() - nanorod.neight3->point().x(),2)+
307			pow(nanorod.neight2->point().y() - nanorod.neight3->point().y(),2)+
308			pow(nanorod.neight2->point().z() - nanorod.neight3->point().z(),2));
309
310			Point_3 point1(nanorod.end1->point().x(), endDist2*endRatio1 + nanorod.neight1->point().y(), nanorod.end1->point().z());
311			Point_3 point2(nanorod.end2->point().x(), endDist2*endRatio2 + nanorod.neight2->point().y(), nanorod.end2->point().z());
312	gezelter	2204	nanorod.end1->point() = point1;
313			nanorod.end2->point() = point2;
314	chuckv	2164
315	gezelter	2204	// Construct normal vector for each face.
316			std::transform( nanoRodPolyhedron.facets_begin(), nanoRodPolyhedron.facets_end(), nanoRodPolyhedron.planes_begin(),
317			Normal_vector());
318	chuckv	2215	}
319			void GeometryBuilder::dumpGeometry(const std::string& geomFileName){
320
321			std::ofstream newGeomFile;
322
323			//create new .md file based on old .md file
324			newGeomFile.open(geomFileName.c_str());
325
326			// Write polyhedron in Object File Format (OFF).
327			CGAL::set_ascii_mode( std::cout);
328			newGeomFile << "OFF" << std::endl << nanoRodPolyhedron.size_of_vertices() << ' '
329			<< nanoRodPolyhedron.size_of_facets() << " 0" << std::endl;
330			std::copy( nanoRodPolyhedron.points_begin(), nanoRodPolyhedron.points_end(),
331			std::ostream_iterator<Point_3>( newGeomFile, "\n"));
332			for ( Facet_iterator i = nanoRodPolyhedron.facets_begin(); i != nanoRodPolyhedron.facets_end(); ++i) {
333			Halfedge_facet_circulator j = i->facet_begin();
334			// Facets in polyhedral surfaces are at least triangles.
335			CGAL_assertion( CGAL::circulator_size(j) >= 3);
336			newGeomFile << CGAL::circulator_size(j) << ' ';
337			do {
338			newGeomFile << ' ' << std::distance(nanoRodPolyhedron.vertices_begin(), j->vertex());
339			} while ( ++j != i->facet_begin());
340			newGeomFile << std::endl;
341			}
342
343			newGeomFile.close();
344
345
346	chuckv	2164
347
348	chuckv	2215	}
349	chuckv	2164
350	chuckv	2215