FEM4EllipticPDE/src/rectilinearstructuredgrid.cpp

// ===========================================================================
// 
//       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();





	}




}