123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433 |
- // ===========================================================================
- //
- // Filename: rectilinearstructuredgrid.cpp
- //
- // Description:
- //
- // Version: 0.0
- // Created: 09/16/2013 01:55:17 PM
- // Revision: none
- // Compiler: Tested with g++
- //
- // Author: M. Andy Kass (MAK)
- //
- // Organisation: US Geological Survey
- //
- //
- // Email: mkass@usgs.gov
- //
- // This program is free software: you can redistribute it and/or modify
- // it under the terms of the GNU General Public License as published by
- // the Free Software Foundation, either version 3 of the License, or
- // (at your option) any later version.
- //
- // This program is distributed in the hope that it will be useful,
- // but WITHOUT ANY WARRANTY; without even the implied warranty of
- // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- // GNU General Public License for more details.
- //
- // You should have received a copy of the GNU General Public License
- // along with this program. If not, see <http://www.gnu.org/licenses/>.
- //
- // ===========================================================================
-
- #include "rectilinearstructuredgrid.h"
-
- namespace formalhaut {
-
- std::ostream &operator<<(std::ostream &stream,
- const RectilinearStructuredGrid &ob) {
-
- return stream;
- }
-
- RectilinearStructuredGrid::RectilinearStructuredGrid() {
- this->TypeOfGrid = rectilinear;
- }
-
- RectilinearStructuredGrid::~RectilinearStructuredGrid() {
-
- }
-
- RectilinearStructuredGrid* RectilinearStructuredGrid::New() {
- return new RectilinearStructuredGrid;
- }
-
- void RectilinearStructuredGrid::Delete() {
- delete this;
- }
-
- void RectilinearStructuredGrid::ComputeCentroids() {
-
- int ntot;
- ntot = this->Nx3(0)*this->Nx3(1)*this->Nx3(2);
- this->VecCentroid.resize(Eigen::NoChange,ntot);
- this->CW.resize(Eigen::NoChange,ntot);
- int mm = 0;
-
- // First z changes, then x then y
- // first depth, then easting, then northing
- //Original
-
- ///*
- //std::cout << this->Dx3v[0].segment(0,1).sum() << std::endl;
- //std::cout << this->Dx3v[0](0) << " " << this->Dx3v[0](1) <<
- // std::endl;
- //std::cout << this->Dx3v[0].segment(0,1).sum() -
- // this->Dx3v[0](1)/2 + this->Dx3v[0](0) +
- // this->Origin(0) << std::endl;
- for (int jj=0;jj<this->Nx3(1);jj++) {
- for (int ii=this->Nx3(0)-1;ii>-1;ii--) {
- //for (int ii=0;ii<this->Nx3(0);ii++) {
- for (int kk=0;kk<this->Nx3(2);kk++) {
-
- //*/
- /*
- for (int ii=0;ii<this->Nx3(0);ii++) {
- for (int jj=0;jj<this->Nx3(1);jj++) {
- for (int kk=0;kk<this->Nx3(2);kk++) {
- */
- this->VecCentroid(0,mm) = this->Origin(0) +
- this->Dx3v[0].segment(0,ii).sum() -
- this->Dx3v[0](ii)/2 + this->Dx3v[0](ii);
- this->VecCentroid(1,mm) = this->Origin(1) +
- this->Dx3v[1].segment(0,jj).sum() -
- this->Dx3v[1](jj)/2 + this->Dx3v[1](jj);
- this->VecCentroid(2,mm) = this->Origin(2) +
- this->Dx3v[2].segment(0,kk).sum() -
- this->Dx3v[2](kk)/2 + this->Dx3v[2](kk);
- this->CW(0,mm) = this->Dx3v[0](ii);
- this->CW(1,mm) = this->Dx3v[1](jj);
- this->CW(2,mm) = this->Dx3v[2](kk);
- mm++;
- }
- }
- }
-
- }
-
- void RectilinearStructuredGrid::SetCellDims(const
- std::vector<VectorXr>& dx3v) {
-
- this->Dx3v[0] = dx3v[0];
- this->Dx3v[1] = dx3v[1];
- this->Dx3v[2] = dx3v[2];
-
- }
-
- std::vector<vtkActor*> RectilinearStructuredGrid::ReturnMeshActor() {
-
- // Coords are edges of cells, not centroids
- vtkFloatArray *xCoords = vtkFloatArray::New();
- vtkFloatArray *yCoords = vtkFloatArray::New();
- vtkFloatArray *zCoords = vtkFloatArray::New();
- Real val;
-
- val = this->Origin(0);
- xCoords->InsertNextValue(val);
- for (int ii=1;ii<=this->Nx3(0)+1;ii++) {
- val = this->Origin(0) + this->Dx3v[0].segment(0,ii-1).sum();
- xCoords->InsertNextValue(val);
- }
- val = this->Origin(1);
- yCoords->InsertNextValue(val);
- for (int ii=1;ii<=this->Nx3(1)+1;ii++) {
- val = this->Origin(1) + this->Dx3v[1].segment(0,ii-1).sum();
- yCoords->InsertNextValue(val);
- }
- val = -1 * this->Origin(2);
- zCoords->InsertNextValue(val);
- for (int ii=1;ii<=this->Nx3(2)+1;ii++) {
- val = this->Origin(2) + this->Dx3v[2].segment(0,ii-1).sum();
- val = -1*val;
- zCoords->InsertNextValue(val);
- }
-
- vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
- rgrid->SetDimensions(this->Nx3(1)+2,this->Nx3(0)+2,this->Nx3(2)+2);
- rgrid->SetXCoordinates(yCoords);
- rgrid->SetYCoordinates(xCoords);
- rgrid->SetZCoordinates(zCoords);
-
- vtkRectilinearGridGeometryFilter *plane =
- vtkRectilinearGridGeometryFilter::New();
- plane->SetInputData(rgrid);
- //plane->SetExtent(0,this->Nx3(0),0,this->Nx3(1),0,this->Nx3(2));
- //Just for testing
- plane->SetExtent(0,0,0,this->Nx3(0)+1,0,this->Nx3(2)+1);
-
- vtkPolyDataMapper *rgridmapper = vtkPolyDataMapper::New();
- rgridmapper->SetInputConnection(plane->GetOutputPort());
-
- vtkActor *wireActor = vtkActor::New();
- wireActor->SetMapper(rgridmapper);
- wireActor->GetProperty()->SetRepresentationToWireframe();
- wireActor->GetProperty()->SetColor(0,1,0);
-
- //Make an actor for each face
- //Second face
- vtkRectilinearGridGeometryFilter *plane1 =
- vtkRectilinearGridGeometryFilter::New();
- plane1->SetInputData(rgrid);
- plane1->SetExtent(this->Nx3(1)+1,this->Nx3(1)+1,0,this->Nx3(0)+1,
- 0,this->Nx3(2)+1);
- vtkPolyDataMapper *rgridmapper1 = vtkPolyDataMapper::New();
- rgridmapper1->SetInputConnection(plane1->GetOutputPort());
- vtkActor *wireActor1 = vtkActor::New();
- wireActor1->SetMapper(rgridmapper1);
- wireActor1->GetProperty()->SetRepresentationToWireframe();
- wireActor1->GetProperty()->SetColor(0,1,0);
- //Third face
- vtkRectilinearGridGeometryFilter *plane2 =
- vtkRectilinearGridGeometryFilter::New();
- plane2->SetInputData(rgrid);
- plane2->SetExtent(0,this->Nx3(1)+1,0,0,
- 0,this->Nx3(2)+1);
- vtkPolyDataMapper *rgridmapper2 = vtkPolyDataMapper::New();
- rgridmapper2->SetInputConnection(plane2->GetOutputPort());
- vtkActor *wireActor2 = vtkActor::New();
- wireActor2->SetMapper(rgridmapper2);
- wireActor2->GetProperty()->SetRepresentationToWireframe();
- wireActor2->GetProperty()->SetColor(0,1,0);
- //Fourth face
- vtkRectilinearGridGeometryFilter *plane3 =
- vtkRectilinearGridGeometryFilter::New();
- plane3->SetInputData(rgrid);
- plane3->SetExtent(0,this->Nx3(1)+1,this->Nx3(0)+1,this->Nx3(0)+1,
- 0,this->Nx3(2)+1);
- vtkPolyDataMapper *rgridmapper3 = vtkPolyDataMapper::New();
- rgridmapper3->SetInputConnection(plane3->GetOutputPort());
- vtkActor *wireActor3 = vtkActor::New();
- wireActor3->SetMapper(rgridmapper3);
- wireActor3->GetProperty()->SetRepresentationToWireframe();
- wireActor3->GetProperty()->SetColor(0,1,0);
- //Fifth face
- vtkRectilinearGridGeometryFilter *plane4 =
- vtkRectilinearGridGeometryFilter::New();
- plane4->SetInputData(rgrid);
- plane4->SetExtent(0,this->Nx3(1)+1,0,this->Nx3(0)+1,
- 0,0);
- vtkPolyDataMapper *rgridmapper4 = vtkPolyDataMapper::New();
- rgridmapper4->SetInputConnection(plane4->GetOutputPort());
- vtkActor *wireActor4 = vtkActor::New();
- wireActor4->SetMapper(rgridmapper4);
- wireActor4->GetProperty()->SetRepresentationToWireframe();
- wireActor4->GetProperty()->SetColor(0,1,0);
- //Sixth face
- vtkRectilinearGridGeometryFilter *plane5 =
- vtkRectilinearGridGeometryFilter::New();
- plane5->SetInputData(rgrid);
- plane5->SetExtent(0,this->Nx3(1)+1,0,this->Nx3(0)+1,
- this->Nx3(2)+1,this->Nx3(2)+1);
- vtkPolyDataMapper *rgridmapper5 = vtkPolyDataMapper::New();
- rgridmapper5->SetInputConnection(plane5->GetOutputPort());
- vtkActor *wireActor5 = vtkActor::New();
- wireActor5->SetMapper(rgridmapper5);
- wireActor5->GetProperty()->SetRepresentationToWireframe();
- wireActor5->GetProperty()->SetColor(0,1,0);
-
- vtkSmartPointer<vtkOutlineFilter> outlineFilter =
- vtkSmartPointer<vtkOutlineFilter>::New();
- outlineFilter->SetInputData(rgrid);
- vtkSmartPointer<vtkDataSetMapper> outlinemapper =
- vtkSmartPointer<vtkDataSetMapper>::New();
- outlinemapper->SetInputConnection(outlineFilter->GetOutputPort());
-
- vtkActor *outlineActor = vtkActor::New();
- outlineActor->SetMapper(outlinemapper);
-
- std::vector<vtkActor*> actorreturn;
- // Order of actors:
- // 1. outlineActor
- // 2-7. plane 0-5
-
- actorreturn.push_back(outlineActor);
- actorreturn.push_back(wireActor);
- actorreturn.push_back(wireActor1);
- actorreturn.push_back(wireActor2);
- actorreturn.push_back(wireActor3);
- actorreturn.push_back(wireActor4);
- actorreturn.push_back(wireActor5);
-
- return actorreturn;
-
-
- }
-
- std::vector<VectorXr> RectilinearStructuredGrid::GetCellDims() {
- return this->Dx3v;
- }
-
- VectorXr RectilinearStructuredGrid::GetCorners(const int& vecpos) {
- VectorXr corners;
- corners.resize(6);
- Vector3r cntr;
-
- cntr = this->VecCentroid.col(vecpos);
-
- corners(0) = cntr(0) - this->CW(0,vecpos)/2;
- corners(1) = cntr(0) + this->CW(0,vecpos)/2;
- corners(2) = cntr(1) - this->CW(1,vecpos)/2;
- corners(3) = cntr(1) + this->CW(1,vecpos)/2;
- corners(4) = cntr(2) - this->CW(2,vecpos)/2;
- corners(5) = cntr(2) + this->CW(2,vecpos)/2;
-
- return corners;
- }
-
- void RectilinearStructuredGrid::DisplayMesh() {
-
- // Coords are edges of cells, not centroids
- vtkFloatArray *xCoords = vtkFloatArray::New();
- vtkFloatArray *yCoords = vtkFloatArray::New();
- vtkFloatArray *zCoords = vtkFloatArray::New();
- Real val;
-
- val = this->Origin(0);
- xCoords->InsertNextValue(val);
- for (int ii=1;ii<=this->Nx3(0)+1;ii++) {
- val = this->Origin(0) + this->Dx3v[0].segment(0,ii-1).sum();
- xCoords->InsertNextValue(val);
- }
- val = this->Origin(1);
- yCoords->InsertNextValue(val);
- for (int ii=1;ii<=this->Nx3(1)+1;ii++) {
- val = this->Origin(1) + this->Dx3v[1].segment(0,ii-1).sum();
- yCoords->InsertNextValue(val);
- }
- val = -1 * this->Origin(2);
- zCoords->InsertNextValue(val);
- for (int ii=1;ii<=this->Nx3(2)+1;ii++) {
- val = this->Origin(2) + this->Dx3v[2].segment(0,ii-1).sum();
- val = -1*val;
- zCoords->InsertNextValue(val);
- }
-
- vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
- rgrid->SetDimensions(this->Nx3(1)+2,this->Nx3(0)+2,this->Nx3(2)+2);
- rgrid->SetXCoordinates(yCoords);
- rgrid->SetYCoordinates(xCoords);
- rgrid->SetZCoordinates(zCoords);
-
- vtkRectilinearGridGeometryFilter *plane =
- vtkRectilinearGridGeometryFilter::New();
- plane->SetInputData(rgrid);
- //plane->SetExtent(0,this->Nx3(0),0,this->Nx3(1),0,this->Nx3(2));
- //Just for testing
- plane->SetExtent(0,0,0,this->Nx3(0)+1,0,this->Nx3(2)+1);
-
- vtkPolyDataMapper *rgridmapper = vtkPolyDataMapper::New();
- rgridmapper->SetInputConnection(plane->GetOutputPort());
-
- vtkActor *wireActor = vtkActor::New();
- wireActor->SetMapper(rgridmapper);
- wireActor->GetProperty()->SetRepresentationToWireframe();
- wireActor->GetProperty()->SetColor(0,1,0);
-
- //Make an actor for each face
- //Second face
- vtkRectilinearGridGeometryFilter *plane1 =
- vtkRectilinearGridGeometryFilter::New();
- plane1->SetInputData(rgrid);
- plane1->SetExtent(this->Nx3(1)+1,this->Nx3(1)+1,0,this->Nx3(0)+1,
- 0,this->Nx3(2)+1);
- vtkPolyDataMapper *rgridmapper1 = vtkPolyDataMapper::New();
- rgridmapper1->SetInputConnection(plane1->GetOutputPort());
- vtkActor *wireActor1 = vtkActor::New();
- wireActor1->SetMapper(rgridmapper1);
- wireActor1->GetProperty()->SetRepresentationToWireframe();
- wireActor1->GetProperty()->SetColor(0,1,0);
- //Third face
- vtkRectilinearGridGeometryFilter *plane2 =
- vtkRectilinearGridGeometryFilter::New();
- plane2->SetInputData(rgrid);
- plane2->SetExtent(0,this->Nx3(1)+1,0,0,
- 0,this->Nx3(2)+1);
- vtkPolyDataMapper *rgridmapper2 = vtkPolyDataMapper::New();
- rgridmapper2->SetInputConnection(plane2->GetOutputPort());
- vtkActor *wireActor2 = vtkActor::New();
- wireActor2->SetMapper(rgridmapper2);
- wireActor2->GetProperty()->SetRepresentationToWireframe();
- wireActor2->GetProperty()->SetColor(0,1,0);
- //Fourth face
- vtkRectilinearGridGeometryFilter *plane3 =
- vtkRectilinearGridGeometryFilter::New();
- plane3->SetInputData(rgrid);
- plane3->SetExtent(0,this->Nx3(1)+1,this->Nx3(0)+1,this->Nx3(0)+1,
- 0,this->Nx3(2)+1);
- vtkPolyDataMapper *rgridmapper3 = vtkPolyDataMapper::New();
- rgridmapper3->SetInputConnection(plane3->GetOutputPort());
- vtkActor *wireActor3 = vtkActor::New();
- wireActor3->SetMapper(rgridmapper3);
- wireActor3->GetProperty()->SetRepresentationToWireframe();
- wireActor3->GetProperty()->SetColor(0,1,0);
- //Fifth face
- vtkRectilinearGridGeometryFilter *plane4 =
- vtkRectilinearGridGeometryFilter::New();
- plane4->SetInputData(rgrid);
- plane4->SetExtent(0,this->Nx3(1)+1,0,this->Nx3(0)+1,
- 0,0);
- vtkPolyDataMapper *rgridmapper4 = vtkPolyDataMapper::New();
- rgridmapper4->SetInputConnection(plane4->GetOutputPort());
- vtkActor *wireActor4 = vtkActor::New();
- wireActor4->SetMapper(rgridmapper4);
- wireActor4->GetProperty()->SetRepresentationToWireframe();
- wireActor4->GetProperty()->SetColor(0,1,0);
- //Sixth face
- vtkRectilinearGridGeometryFilter *plane5 =
- vtkRectilinearGridGeometryFilter::New();
- plane5->SetInputData(rgrid);
- plane5->SetExtent(0,this->Nx3(1)+1,0,this->Nx3(0)+1,
- this->Nx3(2)+1,this->Nx3(2)+1);
- vtkPolyDataMapper *rgridmapper5 = vtkPolyDataMapper::New();
- rgridmapper5->SetInputConnection(plane5->GetOutputPort());
- vtkActor *wireActor5 = vtkActor::New();
- wireActor5->SetMapper(rgridmapper5);
- wireActor5->GetProperty()->SetRepresentationToWireframe();
- wireActor5->GetProperty()->SetColor(0,1,0);
-
- vtkSmartPointer<vtkOutlineFilter> outlineFilter =
- vtkSmartPointer<vtkOutlineFilter>::New();
- outlineFilter->SetInputData(rgrid);
- vtkSmartPointer<vtkDataSetMapper> outlinemapper =
- vtkSmartPointer<vtkDataSetMapper>::New();
- outlinemapper->SetInputConnection(outlineFilter->GetOutputPort());
-
- vtkActor *outlineActor = vtkActor::New();
- outlineActor->SetMapper(outlinemapper);
-
- vtkRenderer *renderer = vtkRenderer::New();
- vtkRenderWindow *renWin = vtkRenderWindow::New();
- renWin->AddRenderer(renderer);
- vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
- iren->SetRenderWindow(renWin);
-
- renderer->AddActor(wireActor);
- //renderer->AddActor(wireActor1);
- renderer->AddActor(wireActor2);
- //renderer->AddActor(wireActor3);
- //renderer->AddActor(wireActor4);
- renderer->AddActor(wireActor5);
- renderer->AddActor(outlineActor);
- renderer->SetBackground(0,0,0);
- renderer->ResetCamera();
- renderer->GetActiveCamera()->Elevation(40.0);
- renderer->GetActiveCamera()->Azimuth(15);
- renderer->GetActiveCamera()->Roll(200);
- renderer->GetActiveCamera()->Zoom(1.0);
-
- renWin->SetSize(900,900);
-
- renWin->Render();
- iren->Start();
-
-
-
-
-
- }
-
-
-
-
- }
|