44 |
|
* |
45 |
|
* Created by Charles F. Vardeman II on 11 Dec 2006. |
46 |
|
* @author Charles F. Vardeman II |
47 |
< |
* @version $Id: ConvexHull.cpp,v 1.10 2008-10-15 18:26:01 chuckv Exp $ |
47 |
> |
* @version $Id: ConvexHull.cpp,v 1.12 2008-10-21 16:44:00 chuckv Exp $ |
48 |
|
* |
49 |
|
*/ |
50 |
|
|
370 |
|
#ifdef HAVE_QHULL |
371 |
|
/* Old options Qt Qu Qg QG0 FA */ |
372 |
|
/* More old opts Qc Qi Pp*/ |
373 |
< |
ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp"), Ns_(200) { |
373 |
> |
ConvexHull::ConvexHull() : Hull(), dim_(3), options_("qhull Qt Pp"), Ns_(200), nTriangles_(0) { |
374 |
|
//If we are doing the mpi version, set up some vectors for data communication |
375 |
|
#ifdef IS_MPI |
376 |
|
|
438 |
|
|
439 |
|
boolT ismalloc = False; |
440 |
|
/* Clean up memory from previous convex hull calculations*/ |
441 |
+ |
|
442 |
|
Triangles_.clear(); |
443 |
|
surfaceSDs_.clear(); |
444 |
|
surfaceSDs_.reserve(Ns_); |
447 |
|
const_cast<char *>(options_.c_str()), NULL, stderr)) { |
448 |
|
|
449 |
|
sprintf(painCave.errMsg, "ConvexHull: Qhull failed to compute convex hull"); |
450 |
< |
painCave.isFatal = 0; |
450 |
> |
painCave.isFatal = 1; |
451 |
|
simError(); |
452 |
|
|
453 |
|
} //qh_new_qhull |
455 |
|
|
456 |
|
#ifdef IS_MPI |
457 |
|
std::vector<double> localPts; |
458 |
+ |
std::vector<double> localVel; |
459 |
|
int localPtArraySize; |
460 |
|
|
461 |
|
|
467 |
|
if (!facet->simplicial){ |
468 |
|
// should never happen with Qt |
469 |
|
sprintf(painCave.errMsg, "ConvexHull: non-simplicaial facet detected"); |
470 |
< |
painCave.isFatal = 0; |
470 |
> |
painCave.isFatal = 1; |
471 |
|
simError(); |
472 |
|
} |
473 |
|
|
519 |
|
localPts.push_back(ptArray[dim_ * idx]); |
520 |
|
localPts.push_back(ptArray[dim_ * idx + 1]); |
521 |
|
localPts.push_back(ptArray[dim_ * idx + 2]); |
522 |
+ |
|
523 |
+ |
Vector3d vel = bodydoubles[idx]->getVel(); |
524 |
+ |
localVel.push_back(vel.x()); |
525 |
+ |
localVel.push_back(vel.y()); |
526 |
+ |
localVel.push_back(vel.z()); |
527 |
+ |
|
528 |
|
localPtsMap.push_back(idx); |
529 |
|
} |
530 |
|
|
543 |
|
int nglobalPts = int(Nsglobal_/3); |
544 |
|
|
545 |
|
|
546 |
< |
std::vector<double> globalPts; |
547 |
< |
globalPts.resize(Nsglobal_); |
546 |
> |
std::vector<double> globalPts(Nsglobal_); |
547 |
> |
std::vector<double> globalVel(Nsglobal_); |
548 |
|
|
549 |
|
isSurfaceID.resize(nglobalPts); |
550 |
|
|
562 |
|
/* gather the potential hull */ |
563 |
|
|
564 |
|
MPI::COMM_WORLD.Allgatherv(&localPts[0],localPtArraySize,MPI::DOUBLE,&globalPts[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE); |
565 |
+ |
MPI::COMM_WORLD.Allgatherv(&localVel[0],localPtArraySize,MPI::DOUBLE,&globalVel[0],&NstoProc_[0],&displs_[0],MPI::DOUBLE); |
566 |
|
|
567 |
|
/* |
568 |
|
if (myrank_ == 0){ |
609 |
|
simError(); |
610 |
|
} //simplicical |
611 |
|
|
612 |
< |
Triangle* face = new Triangle(); |
612 |
> |
Triangle face; |
613 |
|
Vector3d V3dNormal(facet->normal[0],facet->normal[1],facet->normal[2]); |
614 |
< |
face->setNormal(V3dNormal); |
614 |
> |
face.setNormal(V3dNormal); |
615 |
|
|
616 |
|
|
617 |
|
|
618 |
|
RealType faceArea = 0.5*V3dNormal.length(); |
619 |
< |
face->setArea(faceArea); |
619 |
> |
face.setArea(faceArea); |
620 |
|
|
621 |
|
|
622 |
|
vertices = qh_facet3vertex(facet); |
623 |
|
|
624 |
|
coordT *center = qh_getcenter(vertices); |
625 |
|
Vector3d V3dCentroid(center[0], center[1], center[2]); |
626 |
< |
face->setCentroid(V3dCentroid); |
627 |
< |
|
626 |
> |
face.setCentroid(V3dCentroid); |
627 |
> |
Vector3d faceVel = V3Zero; |
628 |
|
FOREACHvertex_(vertices){ |
629 |
|
id = qh_pointid(vertex->point); |
630 |
|
int localindex = id; |
631 |
|
#ifdef IS_MPI |
632 |
< |
|
632 |
> |
Vector3d velVector(globalVel[dim_ * id],globalVel[dim_ * id + 1], globalVel[dim_ * id + 1]); |
633 |
> |
faceVel = faceVel + velVector; |
634 |
|
if (id >= noffset/3 && id < (noffset + localPtArraySize)/3 ){ |
635 |
|
localindex = localPtsMap[id-noffset/3]; |
636 |
+ |
#else |
637 |
+ |
faceVel = faceVel + bodydoubles[localindex]->getVel(); |
638 |
|
#endif |
639 |
< |
face->addVertex(bodydoubles[localindex]); |
639 |
> |
face.addVertex(bodydoubles[localindex]); |
640 |
|
if( !isSurfaceID[id] ){ |
641 |
|
isSurfaceID[id] = true; |
642 |
|
#ifdef IS_MPI |
650 |
|
#ifdef IS_MPI |
651 |
|
|
652 |
|
}else{ |
653 |
< |
face->addVertex(NULL); |
653 |
> |
face.addVertex(NULL); |
654 |
|
} |
655 |
|
#endif |
656 |
|
} //Foreachvertex |
661 |
|
surfaceSDs_.push_back(bodydoubles[id]); |
662 |
|
} |
663 |
|
} |
664 |
< |
|
664 |
> |
*/ |
665 |
> |
face.setFacetVelocity(faceVel/3.0); |
666 |
|
Triangles_.push_back(face); |
667 |
|
qh_settempfree(&vertices); |
668 |
< |
*/ |
668 |
> |
|
669 |
|
} //FORALLfacets |
670 |
|
|
671 |
|
/* |
679 |
|
|
680 |
|
|
681 |
|
Ns_ = surfaceSDs_.size(); |
682 |
+ |
nTriangles_ = Triangles_.size(); |
683 |
|
|
670 |
– |
|
684 |
|
qh_getarea(qh facet_list); |
685 |
|
volume_ = qh totvol; |
686 |
|
area_ = qh totarea; |