// ===========================================================================
//
// 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 .
//
// ===========================================================================
#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;jjNx3(1);jj++) {
for (int ii=this->Nx3(0)-1;ii>-1;ii--) {
//for (int ii=0;iiNx3(0);ii++) {
for (int kk=0;kkNx3(2);kk++) {
//*/
/*
for (int ii=0;iiNx3(0);ii++) {
for (int jj=0;jjNx3(1);jj++) {
for (int kk=0;kkNx3(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& dx3v) {
this->Dx3v[0] = dx3v[0];
this->Dx3v[1] = dx3v[1];
this->Dx3v[2] = dx3v[2];
}
std::vector 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 outlineFilter =
vtkSmartPointer::New();
outlineFilter->SetInputData(rgrid);
vtkSmartPointer outlinemapper =
vtkSmartPointer::New();
outlinemapper->SetInputConnection(outlineFilter->GetOutputPort());
vtkActor *outlineActor = vtkActor::New();
outlineActor->SetMapper(outlinemapper);
std::vector 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 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 outlineFilter =
vtkSmartPointer::New();
outlineFilter->SetInputData(rgrid);
vtkSmartPointer outlinemapper =
vtkSmartPointer::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();
}
}