kernel/Surface.h
Go to the documentation of this file.00001 #ifndef SURFACE_H 00002 #define SURFACE_H 00003 00004 #include <vector> 00005 #include <stdlib.h> 00006 #include <stdio.h> 00007 #include <fstream> 00008 00009 #include "Part.h" 00010 00015 template <class DataType> 00016 class Surface : public Part<DataType,3>, public Facet<DataType> { 00017 public: 00018 typedef Part<DataType,3> PartType; 00019 typedef Connection ConBase; 00020 typedef Segment<DataType,3> SegType; 00021 typedef Facet<DataType> FacetType; 00022 00023 00024 Surface() { 00025 centroid(0) = 0.; centroid(1) = 0.; centroid(2) = 0.; 00026 area = DataType(0.); 00027 volume = DataType(0.); 00028 PartType::DHMat.resize(4,4); 00029 PartType::DHMat.assign(imatrix(4)); 00030 PartType::DHMat_old.resize(4,4); 00031 PartType::DHMat_old.assign(imatrix(4)); 00032 } 00033 ~Surface() {} 00034 00035 virtual void resetPart() { 00036 if (PartType::GetMobility() == 1) { 00037 #ifdef SURFACE_DEBUG_ON 00038 printf("%s->resetPart()\t Surface RESETING cPOINTS\n",PartType::Name().c_str()); 00039 #endif 00040 cPoints.clear(); 00041 std::copy (cPoints_org.begin (), cPoints_org.end (), std::back_inserter (cPoints)); 00042 } 00043 } 00044 00045 virtual DataType Length() const { return 0; } 00046 virtual DataType Area() const { return area; } 00047 virtual DataType Volume() const { return volume; } 00048 virtual void setLength(DataType val) {} 00049 virtual void setArea(DataType val) { area = val; } 00050 virtual void setVolume(DataType val) { volume = val; } 00051 00052 virtual DataType conArea() const { return area; } 00053 void setConArea(DataType val) { area = val; } 00054 00055 virtual bool closure() const { return closed; } 00056 virtual void setClosure(bool val) { closed = val; } 00057 00058 virtual PType Centroid() const { return centroid; } 00059 virtual void setCentroid(PType val) { centroid = val; } 00060 00061 virtual void AddPoint(const PType v) { 00062 cPoints.push_back(v); 00063 } 00064 00065 virtual int GetNumPoints() const { 00066 return ((int) cPoints.size()); 00067 } 00068 00069 virtual PType GetNthPoint(int n) const { 00070 return (cPoints[n]); 00071 } 00072 00073 virtual PType* GetNthPointA(int n) { 00074 return (&cPoints[n]); 00075 } 00076 00077 virtual PType * GetNthPointAddress(int n) { 00078 return &(cPoints[n]); 00079 } 00080 00081 virtual PType * GetNthVelocityAddress(int n) { 00082 return &(velocities_loc[n]); 00083 } 00084 00085 virtual PType UpdateNthPoint(int n, PType val) { 00086 return (cPoints[n] = val); 00087 } 00088 00089 virtual PType GetNthPoint_org(int n) const { 00090 return (cPoints_org[n]); 00091 } 00092 00093 virtual std::vector<PType>& GetCPoints() { 00094 return cPoints; 00095 } 00096 00097 virtual std::vector<PType>& GetCPoints_org() { 00098 return cPoints_org; 00099 } 00100 00101 virtual int GetNumVelocities() const { 00102 return ((int) velocities_loc.size()); 00103 } 00104 00105 virtual PType UpdateNthVelocity(int n, PType val) { 00106 return (cPoints[n] = val); 00107 } 00108 00109 virtual PType GetNthVelocity(int n) const { 00110 return (velocities_loc[n]); 00111 } 00112 00113 virtual std::vector<PType>& GetVelocities() { 00114 return velocities_loc; 00115 } 00116 00117 virtual int GetNumCons() const { 00118 return ((int) connectivity.size()); 00119 } 00120 00121 virtual int GetConnectivity() { return 3;} 00122 virtual void AddConnection(SegType v) {} 00123 virtual void AddConnection(FacetType v) { 00124 connectivity.push_back(v); 00125 } 00126 virtual int GetNumConnections() const { return connectivity.size(); } 00127 virtual ConBase * GetNthConnection(int n) const { return (ConBase*) &connectivity[n]; } 00128 virtual std::vector<FacetType>& GetConnections() { return connectivity; } 00129 virtual const SegType * GetNthSegment(int n) const { return (SegType *) 0; } 00130 virtual const FacetType * GetNthFacet(int n) const { return &connectivity[n]; } 00131 virtual FacetType * EditNthFacet(int n) { return &connectivity[n]; } 00132 virtual SegType * EditNthSegment(int n) { return (SegType *) 0; } 00133 virtual void setNthCon(int i, int val) {} 00134 virtual int getNthCon(int i) const {return -1; } 00135 00136 virtual PType Normal() const { return (PType) 0.; } 00137 virtual void setNormal(PType val) {} 00138 00139 virtual void updatePart(DataType dt,DataType time) { 00140 #ifdef SURFACE_DEBUG_ON 00141 printf("Start Surface update %s dt=%g time=%g\n",PartType::Name().c_str(),dt,time); 00142 #endif 00143 00144 if (PartType::GetMobility() == 1 || time < dt ) { 00145 MType tmp = PartType::GetDHMat(); 00146 #ifdef MOTION_DEBUG_ON 00147 ( mklog() << "\tSurface " << PartType::Name() << " DH_Mat " << tmp << std::endl ).flush(); 00148 ( mklog() << "\n time " << time << " " << PartType::Name() << " DH_MAT\n\t" << (PartType::DHMat-PartType::DHMat_old)*(DataType(1.)/dt) << std::endl ).flush(); 00149 #endif 00150 if ( GetDeformation() == 1 ) { 00151 //deformation 00152 } 00153 prodMVP(tmp,GetCPoints_org(),GetCPoints()); 00154 if (time < dt) 00155 PartType::setDHMat_old(tmp); 00156 numVelMVP(PartType::GetDHMat_old(),PartType::GetDHMat(),GetCPoints(),GetVelocities(),dt); 00157 PartType::setDHMat_old(tmp); 00158 prodMVP(PartType::DHMat,GetCPoints(),GetCPoints()); 00159 measure(); 00160 } 00161 else { 00162 if ( GetDeformation() == 1 ) { 00163 //deformation 00164 } 00165 } 00166 } 00167 00168 virtual void setDeformation(bool val) { PartType::deformable = val; } 00169 virtual bool GetDeformation() const { return PartType::deformable; } 00170 virtual void deformPart(DataType dt,DataType time) {} 00171 00172 virtual void nearestFacet(PType pointNearPart, DataType& minDist, int& facetId, PType& facetC, PType& facetN, int& orientationFlip ) { 00173 Facet<DataType> * facetTmp; 00174 PType a(0.), b(0.); 00175 DataType dist; 00176 00177 for (register int i = 0; i < GetNumConnections(); i++) { 00178 facetTmp = EditNthFacet(i); 00179 facetTmp->setCentroid( DataType(1.)/DataType(3.)*( GetNthPoint(facetTmp->getNthCon(0)) 00180 + GetNthPoint(facetTmp->getNthCon(1)) 00181 + GetNthPoint(facetTmp->getNthCon(2)) ) ); 00182 dist = norm_2( facetTmp->Centroid() - ( pointNearPart ) ); 00183 if (dist < minDist) { 00184 minDist = dist; 00185 facetId = i; 00186 } 00187 } 00188 facetTmp = EditNthFacet(facetId); 00189 facetC = facetTmp->Centroid(); 00190 00191 a = GetNthPoint(facetTmp->getNthCon(1)) - GetNthPoint(facetTmp->getNthCon(0)); 00192 b = GetNthPoint(facetTmp->getNthCon(2)) - GetNthPoint(facetTmp->getNthCon(0)); 00193 facetN = cross3D(a,b); 00194 00195 orientationFlip = 1; 00196 if ( octant(facetN, (pointNearPart - facetC) ) == 0 ) { 00197 facetN = cross3D(b,a); 00198 orientationFlip = 0; 00199 } 00200 00201 } 00202 00203 virtual bool validFacet( const PType& a, const PType& b, const PType& c) { 00204 int valid = 0; 00205 00206 #ifdef SURFACE_DEBUG_ON 00207 ( mklog() << " a " << a << ", b " << b << ", c" << c << std::endl 00208 << " norm_2(a,b)=" << norm_2(a-b) 00209 << " norm_2(a,c)=" << norm_2(a-c) 00210 << " norm_2(b,c)=" << norm_2(b-c) << std::endl ).flush(); 00211 #endif 00212 00213 if (norm_2(a-b) > TOL) { 00214 if (norm_2(a-c) > TOL) { 00215 if (norm_2(b-c) > TOL) { 00216 valid = 1; 00217 } 00218 } 00219 } 00220 #ifdef SURFACE_DEBUG_ON 00221 ( mklog()<< " =Facet " << a << b << c << "valid = " << valid << std::endl ).flush(); 00222 #endif 00223 return ( valid ); 00224 } 00225 00226 virtual void measure() { 00227 area = DataType(0.); 00228 volume = DataType(0.); 00229 centroid(0) = 0.; centroid(1) = 0.; centroid(2) = 0.; 00230 PType atmp(0.); 00231 PType btmp(0.); 00232 PType zero(0.); 00233 int pcheck=1; 00234 00235 #ifdef SURFACE_DEBUG_ON 00236 ( mklog() << "Surface MEASURING " << PartType::Name() << std::endl ).flush(); 00237 #endif 00238 for (register int i = 0; i < GetNumCons(); i++ ) { 00239 atmp = cPoints[connectivity[i].getNthCon(1)]-cPoints[connectivity[i].getNthCon(0)]; 00240 btmp = cPoints[connectivity[i].getNthCon(2)]-cPoints[connectivity[i].getNthCon(0)]; 00241 connectivity[i].setNormal( cross3D( atmp , btmp, pcheck ) ); 00242 00243 connectivity[i].setArea( std::fabs(DataType(0.5)*norm_2( connectivity[i].Normal() )) ); 00244 if (!std::isfinite(connectivity[i].Area())) ( mklog() << "AREA IS A NAN " << i ).flush(); 00245 if (pcheck!=0) { connectivity[i].setNormal(zero); connectivity[i].setArea(0.); } 00246 connectivity[i].setCentroid( DataType(1.)/DataType(3.)*(cPoints[connectivity[i].getNthCon(0)] 00247 + cPoints[connectivity[i].getNthCon(1)] + cPoints[connectivity[i].getNthCon(2)] ) ); 00248 area += connectivity[i].Area(); 00249 volume += innerProd3D(connectivity[i].Centroid(),connectivity[i].Normal()); 00250 00251 #ifdef SURFACE_DEBUG_ON 00252 ( mklog() << "\t\tatmp " << atmp << " btmp " << btmp << " cross3D( atmp , btmp )" 00253 << cross3D( atmp , btmp ) << std::endl ).flush(); 00254 #endif 00255 00256 #ifdef SURFACE_DEBUG_ON 00257 ( mklog() << "Facet centroid " << connectivity[i].Centroid() << " normal " << connectivity[i].Normal() 00258 << " area " << connectivity[i].Area() << " Surface's area " << area << std::endl ).flush(); 00259 #endif 00260 centroid += ( cPoints[connectivity[i].getNthCon(0)] + cPoints[connectivity[i].getNthCon(1)] + cPoints[connectivity[i].getNthCon(2)] ) 00261 *connectivity[i].Area(); 00262 } 00263 centroid *= DataType(1./3.)/area; 00264 if (!std::isfinite(Centroid()(0)) ) ( mklog() << " Centroid0 is NAN " << PartType::Name() ).flush(); 00265 if (!std::isfinite(Centroid()(1)) ) ( mklog() << " Centroid1 is NAN " << PartType::Name() ).flush(); 00266 if (!std::isfinite(Centroid()(2)) ) ( mklog() << " Centroid2 is NAN " << PartType::Name() ).flush(); 00267 volume *= DataType(1./6.); 00268 00269 #ifdef SURFACE_DEBUG_ON 00270 printf("Surface MEASURED %s\n",PartType::Name().c_str()); 00271 if ( PartType::GetMobility() == 1 ) { 00272 ( mklog() << "\tGetMobility() 1 " ).flush(); 00273 } 00274 else { 00275 ( mklog() << "\tGetMobility() 0 " ).flush(); 00276 } 00277 if ( GetDeformation() == 1 ) { 00278 ( mklog() << " GetDeformation() 1\n" ).flush(); 00279 } 00280 else { 00281 ( mklog() << " GetDeformation() 0\n" ).flush(); 00282 00283 } 00284 ( mklog() << "\tPartType::GetCGS() = " << PartType::GetCGS() << std::endl ).flush(); 00285 ( mklog() << "\tarea = " << area << " centroid = " << centroid << std::endl ).flush(); 00286 #endif 00287 00288 } 00289 00290 virtual void read(const int type, std::string file) { 00291 // enum IOType { pointList, VTK_curve, Brep2D, VTK_PL, Brep3D, STL }; 00292 00293 printf("READING %s type = %d dim = %d\n",file.c_str(),type,3); 00294 00295 switch(type) { 00296 case PartType::pointList: { 00297 00298 break; 00299 } 00300 case PartType::VTK_curve: { 00301 00302 break; 00303 } 00304 case PartType::Brep2D: { 00305 00306 break; 00307 } 00308 case PartType::VTK_PL: { 00309 00310 break; 00311 } 00312 case PartType::Brep3D: { 00313 int numPoints, numFacets, dim, rdim, idtmp0, idtmp1, idtmp2; 00314 DataType vtmp; 00315 00316 ( mklog() << "Surface filetype = 3D Brep .iss\n").flush(); 00317 00318 int i=0; 00319 std::string tmp; 00320 00321 std::ifstream infile; 00322 infile.open(file.c_str(), std::ifstream::in); 00323 if (!infile.is_open()) { 00324 ( mklog() <<std::system("pwd")<< "\n\n\n\t\t\tCANNOT OPEN "<<file<<"\n\n\n").flush(); 00325 std::cerr <<std::system("pwd")<< "\n\n\n\t\t\tCANNOT OPEN "<<file<<"\n\n\n"; 00326 exit (-1); 00327 } 00328 00329 PartType::setName(file); 00330 infile >> dim; infile >> rdim; 00331 infile >> numPoints; 00332 while (infile.good() && i < numPoints) { 00333 00334 if (!infile.good()) break; 00335 00336 PType p(0.); 00337 infile >> vtmp; p(0) = vtmp; 00338 infile >> vtmp; p(1) = vtmp; 00339 infile >> vtmp; p(2) = vtmp; 00340 AddPoint(p); 00341 00342 i++; 00343 } 00344 i = 0; 00345 infile >> numFacets; 00346 while (infile.good() && i < numFacets) { 00347 00348 if (!infile.good()) break; 00349 00350 FacetType c; 00351 infile >> idtmp0; c.setNthCon(0,idtmp0); 00352 infile >> idtmp1; c.setNthCon(1,idtmp1); 00353 infile >> idtmp2; c.setNthCon(2,idtmp2); 00354 00355 AddConnection(c); 00356 #ifdef SURFACE_DEBUG_ON 00357 ( mklog()<< PartType::Name()<<" # "<<i<<"\t"<<GetNthFacet(i)->getNthCon(0)<<"\t" 00358 <<GetNthFacet(i)->getNthCon(1)<<"\t"<<GetNthFacet(i)->getNthCon(2)<< std::endl 00359 <<i<<"\t"<<connectivity[i].getNthCon(0)<<"\t"<<connectivity[i].getNthCon(1)<<"\t" 00360 <<connectivity[i].getNthCon(2) ).flush(); 00361 #endif 00362 i++; 00363 } 00364 00365 #ifdef SURFACE_DEBUG_ON 00366 ( mklog() << "GetNthPoint(0) =" << GetNthPoint(0) 00367 << " GetNthPoint(GetNumPoints()-1) =" << GetNthPoint(GetNumPoints()-1) << std::endl 00368 << " ((GetNthPoint(0)-GetNthPoint(GetNumPoints()-1) ) =" 00369 << norm_2(GetNthPoint(0)-GetNthPoint(GetNumPoints()-1) ) << std::endl ).flush(); 00370 #endif 00371 if ( norm_2(GetNthPoint(0)-GetNthPoint(GetNumPoints()-1) ) < TOL ) { 00372 setClosure(1); 00373 } 00374 else 00375 { 00376 setClosure(0); 00377 } 00378 00379 infile.close(); 00380 00381 measure(); 00382 00383 cPoints_org.clear(); 00384 std::copy (cPoints.begin (), cPoints.end (), std::back_inserter (cPoints_org)); 00385 00386 break; 00387 } 00388 case PartType::STL: { 00389 00390 break; 00391 } 00392 } 00393 } 00394 00395 virtual void output(int type, DataType time, FILE * fout) { 00396 ( mklog() << "Writing surface type = "<<type<<" dim ="<<3<<std::endl).flush(); 00397 switch(type) { 00398 case PartType::pointList: { 00399 00400 break; 00401 } 00402 case PartType::VTK_curve: { 00403 00404 break; 00405 } 00406 case PartType::Brep2D: { 00407 00408 break; 00409 } 00410 case PartType::VTK_PL: { 00411 00412 break; 00413 } 00414 case PartType::Brep3D: { 00415 PType a(3); 00416 int totPoints=0; 00417 00418 if (fout == NULL) return; 00419 fprintf(fout,"3 2\n%d\n",GetNumPoints() ); 00420 00421 totPoints = GetNumPoints(); 00422 for (register int j=0; j < totPoints; j++ ) { //GetNumPoints(); j++ ) { 00423 fprintf(fout,"%lf %lf %lf\n",GetNthPoint(j)(0),GetNthPoint(j)(1),GetNthPoint(j)(2) );//cmap[i]->cPoints[j](0),cmap[i]->cPoints[j](1)); 00424 } 00425 00426 //printf("Brep output of %s with numConnections %d\n",PartType::Name().c_str(),GetNumConnections() ); 00427 fprintf(fout,"%d\n",GetNumConnections() ); 00428 for (register int j = 0; j < GetNumConnections(); j++) { 00429 00430 fprintf(fout,"%d\t%d\t%d\n",GetNthFacet(j)->getNthCon(0),GetNthFacet(j)->getNthCon(1),GetNthFacet(j)->getNthCon(2) ); 00431 } 00432 fclose (fout); 00433 00434 break; 00435 } 00436 case PartType::STL: { 00437 if (fout == NULL) return; 00438 fseek (fout, 0, SEEK_END); 00439 00440 fflush(fout); 00441 ( mklog() << "solid "<<PartType::Name()<<std::endl).flush(); 00442 fprintf(fout,"solid %s\n",PartType::Name().c_str()); 00443 for (register int j = 0; j < GetNumConnections(); j++) { 00444 fprintf(fout,"facet normal %g %g %g\nouter loop\n",(GetNthFacet(j)->Normal())(0),(GetNthFacet(j)->Normal())(1),(GetNthFacet(j)->Normal())(2) );//,ni,nj,nk); 00445 if (0) { 00446 fprintf(fout,"vertex %g %g %g\nvertex %g %g %g\nvertex %g %g %g\n", 00447 GetNthPoint(GetNthFacet(j)->getNthCon(0))(2),GetNthPoint(GetNthFacet(j)->getNthCon(0))(1),GetNthPoint(GetNthFacet(j)->getNthCon(0))(0), 00448 GetNthPoint(GetNthFacet(j)->getNthCon(1))(2),GetNthPoint(GetNthFacet(j)->getNthCon(1))(1),GetNthPoint(GetNthFacet(j)->getNthCon(1))(0), 00449 GetNthPoint(GetNthFacet(j)->getNthCon(2))(2),GetNthPoint(GetNthFacet(j)->getNthCon(2))(1),GetNthPoint(GetNthFacet(j)->getNthCon(2))(0) ); 00450 00451 } 00452 else { 00453 fprintf(fout,"vertex %g %g %g\nvertex %g %g %g\nvertex %g %g %g\n", 00454 GetNthPoint(GetNthFacet(j)->getNthCon(0))(0),GetNthPoint(GetNthFacet(j)->getNthCon(0))(1),GetNthPoint(GetNthFacet(j)->getNthCon(0))(2), 00455 GetNthPoint(GetNthFacet(j)->getNthCon(1))(0),GetNthPoint(GetNthFacet(j)->getNthCon(1))(1),GetNthPoint(GetNthFacet(j)->getNthCon(1))(2), 00456 GetNthPoint(GetNthFacet(j)->getNthCon(2))(0),GetNthPoint(GetNthFacet(j)->getNthCon(2))(1),GetNthPoint(GetNthFacet(j)->getNthCon(2))(2) ); 00457 } 00458 fprintf(fout,"endloop\nendfacet\n"); 00459 } 00460 fprintf(fout,"endsolid %s\n",PartType::Name().c_str()); 00461 00462 break; 00463 } 00464 case PartType::VelDat: { 00465 if (fout == NULL) return; 00466 00467 fflush(fout); 00468 ( mklog() << "Writing velDat " << " GetNumPoints() = " << GetNumPoints() << " GetNumVelocities() = " << GetNumVelocities() << std::endl ).flush(); 00469 for (register int j = 0; j < GetNumPoints(); j++) { 00470 fprintf(fout,"%g %g %g\n",(GetNthVelocity(j))(0),(GetNthVelocity(j))(1),(GetNthVelocity(j))(2) ); 00471 } 00472 fclose (fout); 00473 break; 00474 } 00475 } 00476 } 00477 00478 virtual void output(int type, DataType time, std::string file) { 00479 ( mklog() << "Writing surface type = "<<type<<" dim ="<<3<<std::endl).flush(); 00480 switch(type) { 00481 case PartType::pointList: { 00482 00483 break; 00484 } 00485 case PartType::VTK_curve: { 00486 00487 break; 00488 } 00489 case PartType::Brep2D: { 00490 00491 break; 00492 } 00493 case PartType::VTK_PL: { 00494 00495 break; 00496 } 00497 case PartType::Brep3D: { 00498 PType a(3); 00499 int totPoints=0; 00500 00501 FILE * fout; 00502 fout = fopen (file.c_str(),"w"); 00503 if (fout == NULL) return; 00504 fprintf(fout,"3 2\n%d\n",GetNumPoints() ); 00505 00506 totPoints = GetNumPoints(); 00507 for (register int j=0; j < totPoints; j++ ) { //GetNumPoints(); j++ ) { 00508 fprintf(fout,"%lf %lf %lf\n",GetNthPoint(j)(0),GetNthPoint(j)(1),GetNthPoint(j)(2) );//cmap[i]->cPoints[j](0),cmap[i]->cPoints[j](1)); 00509 } 00510 00511 printf("Brep output of %s with numConnections %d\n",PartType::Name().c_str(),GetNumConnections() ); 00512 fprintf(fout,"%d\n",GetNumConnections() ); 00513 for (register int j = 0; j < GetNumConnections(); j++) { 00514 00515 fprintf(fout,"%d\t%d\t%d\n",GetNthFacet(j)->getNthCon(0),GetNthFacet(j)->getNthCon(1),GetNthFacet(j)->getNthCon(2) ); 00516 } 00517 fclose (fout); 00518 00519 break; 00520 } 00521 case PartType::STL: { 00522 FILE * fout; 00523 fout = fopen (file.c_str(),"a+"); 00524 if (fout == NULL) return; 00525 00526 fflush(fout); 00527 fprintf(fout,"solid %s\n",PartType::Name().c_str()); 00528 for (register int j = 0; j < GetNumConnections(); j++) { 00529 fprintf(fout,"facet normal %g %g %g\nouter loop\n",(GetNthFacet(j)->Normal())(0),(GetNthFacet(j)->Normal())(1),(GetNthFacet(j)->Normal())(2) );//,ni,nj,nk); 00530 if (0) { 00531 fprintf(fout,"vertex %g %g %g\nvertex %g %g %g\nvertex %g %g %g\n", 00532 GetNthPoint(GetNthFacet(j)->getNthCon(0))(2),GetNthPoint(GetNthFacet(j)->getNthCon(0))(1),GetNthPoint(GetNthFacet(j)->getNthCon(0))(0), 00533 GetNthPoint(GetNthFacet(j)->getNthCon(1))(2),GetNthPoint(GetNthFacet(j)->getNthCon(1))(1),GetNthPoint(GetNthFacet(j)->getNthCon(1))(0), 00534 GetNthPoint(GetNthFacet(j)->getNthCon(2))(2),GetNthPoint(GetNthFacet(j)->getNthCon(2))(1),GetNthPoint(GetNthFacet(j)->getNthCon(2))(0) ); 00535 00536 } 00537 else { 00538 fprintf(fout,"vertex %g %g %g\nvertex %g %g %g\nvertex %g %g %g\n", 00539 GetNthPoint(GetNthFacet(j)->getNthCon(0))(0),GetNthPoint(GetNthFacet(j)->getNthCon(0))(1),GetNthPoint(GetNthFacet(j)->getNthCon(0))(2), 00540 GetNthPoint(GetNthFacet(j)->getNthCon(1))(0),GetNthPoint(GetNthFacet(j)->getNthCon(1))(1),GetNthPoint(GetNthFacet(j)->getNthCon(1))(2), 00541 GetNthPoint(GetNthFacet(j)->getNthCon(2))(0),GetNthPoint(GetNthFacet(j)->getNthCon(2))(1),GetNthPoint(GetNthFacet(j)->getNthCon(2))(2) ); 00542 } 00543 fprintf(fout,"endloop\nendfacet\n"); 00544 } 00545 fprintf(fout,"endsolid %s\n",PartType::Name().c_str()); 00546 fclose (fout); 00547 00548 break; 00549 } 00550 case PartType::VelDat: { 00551 FILE * fout; 00552 fout = fopen (file.c_str(),"a+"); 00553 if (fout == NULL) return; 00554 00555 fflush(fout); 00556 ( mklog() << "Writing velDat " << file << " GetNumPoints() = " << GetNumPoints() << " GetNumVelocities() = " << GetNumVelocities() << std::endl ).flush(); 00557 for (register int j = 0; j < GetNumPoints(); j++) { 00558 fprintf(fout,"%g %g %g\n",(GetNthVelocity(j))(0),(GetNthVelocity(j))(1),(GetNthVelocity(j))(2) ); 00559 00560 } 00561 fclose (fout); 00562 00563 break; 00564 } 00565 } 00566 } 00567 00568 virtual void Restart(std::ifstream& ifs, int& pos, double& t, double& dt) { 00569 #ifdef CHECKRESTART_DEBUG_ON 00570 ( mklog() << "Restarting Surface " << PartType::Name() << " pos " << pos << " t = " << t << " dt = " << dt << std::endl ).flush(); 00571 #endif 00572 ifs.seekg(pos); 00573 ifs.read((char*)¢roid(0),sizeof(DataType)); 00574 ifs.read((char*)¢roid(1),sizeof(DataType)); 00575 ifs.read((char*)¢roid(2),sizeof(DataType)); 00576 ifs.read((char*)&area,sizeof(DataType)); 00577 ifs.read((char*)&volume,sizeof(DataType)); 00578 00579 ifs.read((char*)&num_cPoints_org,sizeof(int)); 00580 ifs.read((char*)&num_cPoints,sizeof(int)); 00581 ifs.read((char*)&num_connections,sizeof(int)); 00582 ifs.read((char*)&num_velocities,sizeof(int)); 00583 00584 for (register int i=0; i < num_cPoints_org; i++) { 00585 ifs.read((char*)&cPoints_org[i](0),sizeof(DataType)); 00586 ifs.read((char*)&cPoints_org[i](1),sizeof(DataType)); 00587 ifs.read((char*)&cPoints_org[i](2),sizeof(DataType)); 00588 } 00589 for (register int i=0; i < num_cPoints; i++) { 00590 ifs.read((char*)&cPoints[i](0),sizeof(DataType)); 00591 ifs.read((char*)&cPoints[i](1),sizeof(DataType)); 00592 ifs.read((char*)&cPoints[i](2),sizeof(DataType)); 00593 } 00594 for (register int i=0; i < num_connections; i++) { 00595 ifs.read((char*)&(connectivity[i].cPid[0]),sizeof(int)); 00596 ifs.read((char*)&(connectivity[i].cPid[1]),sizeof(int)); 00597 ifs.read((char*)&(connectivity[i].cPid[2]),sizeof(int)); 00598 } 00599 for (register int i=0; i < num_velocities; i++) { 00600 ifs.read((char*)&velocities_loc[i](0),sizeof(DataType)); 00601 ifs.read((char*)&velocities_loc[i](1),sizeof(DataType)); 00602 ifs.read((char*)&velocities_loc[i](2),sizeof(DataType)); 00603 } 00604 00605 ifs.read((char*)&netP(0),sizeof(DataType)); 00606 ifs.read((char*)&netP(1),sizeof(DataType)); 00607 ifs.read((char*)&netP(2),sizeof(DataType)); 00608 ifs.read((char*)&momP(0),sizeof(DataType)); 00609 ifs.read((char*)&momP(1),sizeof(DataType)); 00610 ifs.read((char*)&momP(2),sizeof(DataType)); 00611 ifs.read((char*)&xAxis(0),sizeof(DataType)); 00612 ifs.read((char*)&xAxis(1),sizeof(DataType)); 00613 ifs.read((char*)&xAxis(2),sizeof(DataType)); 00614 ifs.read((char*)&yAxis(0),sizeof(DataType)); 00615 ifs.read((char*)&yAxis(1),sizeof(DataType)); 00616 ifs.read((char*)&yAxis(2),sizeof(DataType)); 00617 ifs.read((char*)&zAxis(0),sizeof(DataType)); 00618 ifs.read((char*)&zAxis(1),sizeof(DataType)); 00619 ifs.read((char*)&zAxis(2),sizeof(DataType)); 00620 00621 pos = ifs.tellg(); 00622 00623 PartType::RestartPart(ifs,pos,t,dt); 00624 pos = ifs.tellg(); 00625 } 00626 00627 virtual void Checkpointing(std::ofstream& ofs) { 00628 #ifdef CHECKRESTART_DEBUG_ON 00629 ( mklog() << "\nCHECKPOINTING Surface " << PartType::Name() << std::endl ).flush(); 00630 #endif 00631 ofs.write((char*)¢roid(0),sizeof(DataType)); 00632 ofs.write((char*)¢roid(1),sizeof(DataType)); 00633 ofs.write((char*)¢roid(2),sizeof(DataType)); 00634 ofs.write((char*)&area,sizeof(DataType)); 00635 ofs.write((char*)&volume,sizeof(DataType)); 00636 00637 num_cPoints_org=cPoints_org.size(); 00638 num_cPoints=cPoints.size(); 00639 num_connections=connectivity.size(); 00640 num_velocities=velocities_loc.size(); 00641 ofs.write((char*)&num_cPoints_org,sizeof(int)); 00642 ofs.write((char*)&num_cPoints,sizeof(int)); 00643 ofs.write((char*)&num_connections,sizeof(int)); 00644 ofs.write((char*)&num_velocities,sizeof(int)); 00645 for (register int i=0; i < num_cPoints_org; i++) { 00646 ofs.write((char*)&cPoints_org[i](0),sizeof(DataType)); 00647 ofs.write((char*)&cPoints_org[i](1),sizeof(DataType)); 00648 ofs.write((char*)&cPoints_org[i](2),sizeof(DataType)); 00649 } 00650 for (register int i=0; i < num_cPoints; i++) { 00651 ofs.write((char*)&cPoints[i](0),sizeof(DataType)); 00652 ofs.write((char*)&cPoints[i](1),sizeof(DataType)); 00653 ofs.write((char*)&cPoints[i](2),sizeof(DataType)); 00654 } 00655 for (register int i=0; i < num_connections; i++) { 00656 ofs.write((char*)&(connectivity[i].cPid[0]),3*sizeof(int)); 00657 } 00658 for (register int i=0; i < num_velocities; i++) { 00659 ofs.write((char*)&velocities_loc[i](0),sizeof(DataType)); 00660 ofs.write((char*)&velocities_loc[i](1),sizeof(DataType)); 00661 ofs.write((char*)&velocities_loc[i](2),sizeof(DataType)); 00662 } 00663 00664 ofs.write((char*)&netP(0),sizeof(DataType)); 00665 ofs.write((char*)&netP(1),sizeof(DataType)); 00666 ofs.write((char*)&netP(2),sizeof(DataType)); 00667 ofs.write((char*)&momP(0),sizeof(DataType)); 00668 ofs.write((char*)&momP(1),sizeof(DataType)); 00669 ofs.write((char*)&momP(2),sizeof(DataType)); 00670 ofs.write((char*)&xAxis(0),sizeof(DataType)); 00671 ofs.write((char*)&xAxis(1),sizeof(DataType)); 00672 ofs.write((char*)&xAxis(2),sizeof(DataType)); 00673 ofs.write((char*)&yAxis(0),sizeof(DataType)); 00674 ofs.write((char*)&yAxis(1),sizeof(DataType)); 00675 ofs.write((char*)&yAxis(2),sizeof(DataType)); 00676 ofs.write((char*)&zAxis(0),sizeof(DataType)); 00677 ofs.write((char*)&zAxis(1),sizeof(DataType)); 00678 ofs.write((char*)&zAxis(2),sizeof(DataType)); 00679 00680 PartType::CheckpointPart(ofs); 00681 } 00682 00683 virtual void updateAxis() { 00684 PType atmp; 00685 MType mtmp(PartType::GetDHMat()); 00686 mtmp(0,3) = 0; mtmp(1,3) = 0; mtmp(2,3) = 0; 00687 atmp(0) = 1.; atmp(1) = 0.; atmp(2) = 0.; 00688 xAxis = prodMP(mtmp,atmp); 00689 xAxis = normalize(xAxis); 00690 atmp(0) = 0.; atmp(1) = 1.; atmp(2) = 0.; 00691 yAxis = prodMP(mtmp,atmp); 00692 yAxis = normalize(yAxis); 00693 atmp(0) = 0.; atmp(1) = 0.; atmp(2) = 1.; 00694 zAxis = prodMP(mtmp,atmp); 00695 zAxis = normalize(zAxis); 00696 } 00697 00698 virtual void pressureForce( multi_index_type* cons, DataType* press, DataType scale ) { 00699 DataType ptmp = 0., stmp = 0.; 00700 PType rtmp; 00701 int offset = PartType::GetCGS(); 00702 00703 num_connections = connectivity.size(); 00704 pressureForce_con.resize(num_connections); 00705 netP = 0.; 00706 momP = 0.; 00707 stmp = DataType(0.5)*scale*scale; 00708 00709 updateAxis(); 00710 00711 for (register int i = 0; i < num_connections; i++) { 00712 ptmp = DataType(1.)/DataType(3.) * ( press[cons[offset+i](0)] + press[cons[offset+i](1)] + press[cons[offset+i](2)] ); 00713 pressureForce_con[i] = DataType(-1.) * connectivity[i].Normal() * ptmp * stmp; 00714 netP += pressureForce_con[i]; 00715 00716 rtmp = connectivity[i].Centroid() - Centroid(); 00717 momP += cross3D(rtmp,pressureForce_con[i]); 00718 } 00719 } 00720 00721 virtual void logLoad( std::ofstream& ofs, DataType time_, int step ) { 00722 updateAxis(); 00723 ofs << step << "\t" << time_ << "\t" << PartType::Name() << "\t" 00724 << centroid << "\t" << xAxis << "\t" << yAxis << "\t" << zAxis << "\t" << netP << "\t" << momP << "\t" 00725 << projectScalar(netP,xAxis) << "\t" << projectScalar(netP,yAxis) << "\t" << projectScalar(netP,zAxis) << "\t" //<< std::endl; 00726 << "\t" << area << "\t" << volume << "\n"; 00727 } 00728 00729 virtual PType GetNetP() { return netP; } 00730 virtual PType GetMomP() { return momP; } 00731 00732 virtual PType GetXaxis() { return xAxis; } 00733 virtual PType GetYaxis() { return yAxis; } 00734 virtual PType GetZaxis() { return zAxis; } 00735 00736 protected: 00737 bool closed; 00738 bool loft; 00739 PType centroid; 00740 PType netP, momP, xAxis, yAxis, zAxis; 00741 DataType area; 00742 DataType volume; 00743 std::vector<PType> cPoints; 00744 std::vector<PType> cPoints_org; 00745 std::vector<PType> velocities_loc; 00746 std::vector<PType> pressureForce_con; 00747 std::vector<FacetType> connectivity; 00748 int num_cPoints_org, num_cPoints, num_connections, num_velocities; 00749 }; 00750 00751 #endif // SURFACE_H
Generated on Thu Jun 30 02:18:58 2016 for Motion - rigid body dynamics solver by
1.6.3