Browse Source

Initial creation of meshtools reader for FEM4EllipticPDE. Corrected the initial name and email.

readMeshtools
M. Andy Kass 8 years ago
parent
commit
1a7ac7bfeb

+ 7
- 0
examples/utReadRectilinearMesh.cpp View File

1
+#include "LemmaObject.h"
2
+
3
+int main(int argc, char**argv) {
4
+
5
+    std::cout << "Hello World!" << std::endl;
6
+
7
+}

+ 47
- 0
include/exceptions.h View File

1
+// ===========================================================================
2
+// 
3
+//       Filename:  exceptions.h
4
+// 
5
+//    Description:  Exception classes
6
+// 
7
+//        Version:  0.0
8
+//        Created:  03/19/2013 12:22:19 PM
9
+//       Revision:  none
10
+//       Compiler:  Tested with g++ 
11
+// 
12
+//         Author:  M. Andy Kass (MAK)
13
+//
14
+//   Organisation:  US Geological Survey
15
+//                  
16
+//
17
+//          Email:  mkass@usgs.gov
18
+// 
19
+//  This program is free software: you can redistribute it and/or modify
20
+//  it under the terms of the GNU General Public License as published by
21
+//  the Free Software Foundation, either version 3 of the License, or
22
+//  (at your option) any later version.
23
+//
24
+//  This program is distributed in the hope that it will be useful,
25
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
26
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27
+//  GNU General Public License for more details.
28
+//
29
+//  You should have received a copy of the GNU General Public License
30
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
31
+//
32
+// ===========================================================================
33
+
34
+#ifndef __exceptions_h
35
+#define __exceptions_h
36
+
37
+#include "formalhaut.h"
38
+
39
+namespace formalhaut {
40
+
41
+/// Generic file I/O error
42
+class GenericFileIOError : public std::runtime_error {
43
+	public: GenericFileIOError(const std::string &filename);
44
+};
45
+
46
+} //end of namespace formalhaut
47
+#endif

+ 87
- 0
include/grid.h View File

1
+// ===========================================================================
2
+//
3
+//       Filename:  grid.h
4
+//
5
+//    Description:
6
+//
7
+//        Version:  0.0
8
+//        Created:  10/25/2009 07:38:32 AM
9
+//       Revision:  none
10
+//       Compiler:  g++ (c++)
11
+//
12
+//         Author:  Trevor Irons (ti)
13
+//         			M. Andy Kass (MAK)
14
+//
15
+//   Organisation:  Colorado School of Mines (CSM)
16
+//                  United States Geological Survey (USGS)
17
+//
18
+//          Email:  tirons@mines.edu, tirons@usgs.gov
19
+//
20
+//  This program is free software: you can redistribute it and/or modify
21
+//  it under the terms of the GNU General Public License as published by
22
+//  the Free Software Foundation, either version 3 of the License, or
23
+//  (at your option) any later version.
24
+//
25
+//  This program is distributed in the hope that it will be useful,
26
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
27
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
28
+//  GNU General Public License for more details.
29
+//
30
+//  You should have received a copy of the GNU General Public License
31
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
32
+//
33
+// ===========================================================================
34
+
35
+#ifndef __GRID_H
36
+#define __GRID_H
37
+#include "formalhaut.h"
38
+
39
+namespace formalhaut {
40
+
41
+// ===============================================================================
42
+//        Class:  Grid
43
+/// \brief Abstract base class for grids.
44
+/// \details Grids define values in space, they may either be an integral part
45
+/// of models, or simply the output of a forward modelling procedure.
46
+// ===============================================================================
47
+
48
+class Grid {
49
+
50
+    friend std::ostream &operator<<(std::ostream &stream,
51
+            const Grid &ob);
52
+
53
+    public:
54
+
55
+    // ====================  LIFECYCLE     ===================================
56
+
57
+    // ====================  OPERATORS     ===================================
58
+
59
+    // ====================  OPERATIONS    ===================================
60
+
61
+    // ====================  ACCESS        ===================================
62
+
63
+    // ====================  INQUIRY       ===================================
64
+
65
+	gridtype GetGridType();
66
+
67
+    protected:
68
+
69
+    // ====================  LIFECYCLE     ===================================
70
+
71
+    /// Default protected constructor.
72
+    Grid ();
73
+
74
+    /// Default protected constructor.
75
+    ~Grid ();
76
+
77
+    // ====================  DATA MEMBERS  ===================================
78
+
79
+	gridtype TypeOfGrid;
80
+
81
+    private:
82
+
83
+}; // -----  end of class  Grid  -----
84
+
85
+
86
+} // namespace formalhaut
87
+#endif

+ 117
- 0
include/logging.h View File

1
+// ===========================================================================
2
+// 
3
+//       Filename:  logging.h
4
+// 
5
+//    Description:  Simple logging class.
6
+// 
7
+//        Version:  0.0
8
+//        Created:  09/12/2013 03:23:30 PM
9
+//       Revision:  none
10
+//       Compiler:  Tested with g++ 
11
+// 
12
+//         Author:  M. Andy Kass (MAK)
13
+//
14
+//   Organisation:  US Geological Survey
15
+//                  
16
+//
17
+//          Email:  mkass@usgs.gov
18
+// 
19
+//  This program is free software: you can redistribute it and/or modify
20
+//  it under the terms of the GNU General Public License as published by
21
+//  the Free Software Foundation, either version 3 of the License, or
22
+//  (at your option) any later version.
23
+//
24
+//  This program is distributed in the hope that it will be useful,
25
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
26
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27
+//  GNU General Public License for more details.
28
+//
29
+//  You should have received a copy of the GNU General Public License
30
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
31
+//
32
+// ===========================================================================
33
+
34
+#ifndef __logging_h
35
+#define __logging_h
36
+
37
+#include "formalhaut.h"
38
+#include "exceptions.h"
39
+#include <time.h>
40
+
41
+namespace formalhaut {
42
+
43
+
44
+// ===================================================================
45
+//        Class:  Logging
46
+/// \brief  
47
+/// \details 
48
+// ===================================================================
49
+class Logging {
50
+
51
+	friend std::ostream &operator<<(std::ostream &stream,
52
+		const Logging &ob);
53
+
54
+public:
55
+
56
+	// ====================  LIFECYCLE     =======================
57
+
58
+	static Logging* New();                    
59
+
60
+	void Delete();
61
+
62
+	// ====================  OPERATORS     =======================
63
+
64
+	// ====================  OPERATIONS    =======================
65
+
66
+	/// Close the log file
67
+	void CloseLogFile();
68
+
69
+	/// Write a single line to the log file
70
+	void WriteLogLine(const std::string line);
71
+
72
+	/// Write header information in the logger
73
+	void WriteHeader(const std::string& name,const std::string& ver);
74
+
75
+	/// Get current date and time, formatted
76
+	std::string GetDateTime();
77
+
78
+	/// Write successful completion
79
+	void WriteCloser();
80
+
81
+	/// Write a long dashed line
82
+	void WriteDashedLine();
83
+
84
+	// ====================  ACCESS        =======================
85
+
86
+	/// Set and open the log file
87
+	void SetLogFile(const std::string& fname);
88
+
89
+	// ====================  INQUIRY       =======================
90
+
91
+	std::string GetLogFile();
92
+
93
+
94
+protected:
95
+
96
+	// ====================  LIFECYCLE     =======================
97
+
98
+	/// Default protected constructor. 
99
+	Logging ();                    
100
+
101
+	/// Default protected constructor. 
102
+	~Logging ();                    
103
+
104
+	// ====================  DATA MEMBERS  =========================
105
+
106
+	// Log file name
107
+	std::string LogFile;
108
+
109
+	// File object
110
+	std::ofstream outfile;
111
+
112
+private:
113
+
114
+}; // -----  end of class  Logging  -----
115
+
116
+} // end of namespace
117
+#endif

+ 110
- 0
include/rectilinearstructuredgrid.h View File

1
+// ===========================================================================
2
+// 
3
+//       Filename:  rectilinearstructuredgrid.h
4
+// 
5
+//    Description:  
6
+// 
7
+//        Version:  0.0
8
+//        Created:  09/16/2013 01:49:40 PM
9
+//       Revision:  none
10
+//       Compiler:  Tested with g++ 
11
+// 
12
+//         Author:  M. Andy Kass (MAK)
13
+//
14
+//   Organisation:  US Geological Survey
15
+//                  
16
+//
17
+//          Email:  mkass@usgs.gov
18
+// 
19
+//  This program is free software: you can redistribute it and/or modify
20
+//  it under the terms of the GNU General Public License as published by
21
+//  the Free Software Foundation, either version 3 of the License, or
22
+//  (at your option) any later version.
23
+//
24
+//  This program is distributed in the hope that it will be useful,
25
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
26
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27
+//  GNU General Public License for more details.
28
+//
29
+//  You should have received a copy of the GNU General Public License
30
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
31
+//
32
+// ===========================================================================
33
+
34
+#ifndef __RECTILINEARSTRUCTUREDGRID_H
35
+#define __RECTILINEARSTRUCTUREDGRID_H
36
+#include "formalhaut.h"
37
+#include "structuredgrid.h"
38
+#include "exceptions.h"
39
+
40
+namespace formalhaut {
41
+
42
+
43
+// ===================================================================
44
+//        Class:  RectilinearStructuredGrid
45
+/// \brief  
46
+/// \details 
47
+// ===================================================================
48
+class RectilinearStructuredGrid : public StructuredGrid {
49
+	friend std::ostream &operator<<(std::ostream &stream,
50
+		const RectilinearStructuredGrid &ob);
51
+
52
+public:
53
+
54
+	// ====================  LIFECYCLE     =======================
55
+
56
+	static RectilinearStructuredGrid* New();                    
57
+
58
+	void Delete();
59
+
60
+	// ====================  OPERATORS     =======================
61
+
62
+	// ====================  OPERATIONS    =======================
63
+	
64
+	/// Compute the centroid of each vectorized cell
65
+	void ComputeCentroids();
66
+
67
+	/// Set the dimensions of the cells
68
+	void SetCellDims(const std::vector<VectorXr>& dx3v);
69
+
70
+	// ====================  ACCESS        =======================
71
+	
72
+	/// Return actors containing the mesh (7 actors returned.  1st is 
73
+	/// bounding box and the next 6 are prism faces)
74
+	std::vector<vtkActor*> ReturnMeshActor();
75
+
76
+	// ====================  INQUIRY       =======================
77
+	
78
+	/// Return dx,dy,dz
79
+	std::vector<VectorXr> GetCellDims();
80
+
81
+	/// Return a particular set of corners
82
+	VectorXr GetCorners(const int& vecpos);
83
+
84
+	/// Display the grid
85
+	void DisplayMesh();
86
+
87
+
88
+
89
+protected:
90
+
91
+	// ====================  LIFECYCLE     =======================
92
+
93
+	/// Default protected constructor. 
94
+	RectilinearStructuredGrid ();                    
95
+
96
+	/// Default protected constructor. 
97
+	~RectilinearStructuredGrid ();  
98
+
99
+	/// Vectorized cell widths.  Not the most memory efficient, but 
100
+	/// efficient computationally (by that I mean I'm lazy and tired
101
+	/// and this seems like the easiest way).
102
+	Vector3Xr CW;	
103
+
104
+	// ====================  DATA MEMBERS  =========================
105
+
106
+private:
107
+
108
+}; // -----  end of class  RectilinearStructuredGrid  -----
109
+} // end of namespace formalhaut
110
+#endif

+ 144
- 0
include/structuredgrid.h View File

1
+// ===========================================================================
2
+// 
3
+//       Filename:  structuredgrid.h
4
+// 
5
+//    Description:  Derived class for structured grids.  Base class to 
6
+//    				regular and rectilinear grids.
7
+// 
8
+//        Version:  0.0
9
+//        Created:  09/12/2013 08:29:26 AM
10
+//       Revision:  none
11
+//       Compiler:  Tested with g++ 
12
+// 
13
+//         Author:  M. Andy Kass (MAK)
14
+//
15
+//   Organisation:  US Geological Survey
16
+//                  
17
+//
18
+//          Email:  mkass@usgs.gov
19
+// 
20
+//  This program is free software: you can redistribute it and/or modify
21
+//  it under the terms of the GNU General Public License as published by
22
+//  the Free Software Foundation, either version 3 of the License, or
23
+//  (at your option) any later version.
24
+//
25
+//  This program is distributed in the hope that it will be useful,
26
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
27
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
28
+//  GNU General Public License for more details.
29
+//
30
+//  You should have received a copy of the GNU General Public License
31
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
32
+//
33
+// ===========================================================================
34
+
35
+#ifndef __STRUCTUREDGRID_H
36
+#define __STRUCTUREDGRID_H
37
+#include "formalhaut.h"
38
+#include "grid.h"
39
+#include "exceptions.h"
40
+#include "logging.h"
41
+
42
+namespace formalhaut {
43
+
44
+// ===================================================================
45
+//        Class:  StructuredGrid
46
+/// \brief  
47
+/// \details 
48
+// ===================================================================
49
+class StructuredGrid : public Grid {
50
+	friend std::ostream &operator<<(std::ostream &stream,
51
+		const StructuredGrid &ob);
52
+
53
+public:
54
+
55
+	// ====================  LIFECYCLE     =======================
56
+
57
+	static StructuredGrid* New();                    
58
+
59
+	void Delete();
60
+
61
+	// ====================  OPERATORS     =======================
62
+
63
+	// ====================  OPERATIONS    =======================
64
+
65
+	/// Set the number of cells in each direction
66
+	void SetNCells(const Eigen::Matrix<int,3,1>& ncells);
67
+
68
+	/// Set origin (top south west)
69
+	void SetOrigin(const Vector3r& origin);
70
+
71
+
72
+	// ====================  ACCESS        =======================
73
+
74
+	/// Read mesh parameters from Meshtools format
75
+	void ReadMeshtoolsGrid(const std::string& fname);
76
+
77
+	/// Attach a log object for writing lines to log file
78
+	void AttachLogObject(Logging* logobj);
79
+
80
+	// ====================  INQUIRY       =======================
81
+
82
+	/// Given a cell count in each direction, return the position that
83
+	/// cell has in the unwrapped vector.  Indexed from zero.
84
+	int GetVectorPosition(const Eigen::Matrix<int,3,1>& cellpos);
85
+
86
+	/// Given a position in space, determine if the position exists
87
+	/// within the model region and then if so, in what cell does that 
88
+	/// position exist.
89
+	int QueryCellPosition(const Vector3r& cellpos);
90
+
91
+	/// Return vector with the number of cells in each direction
92
+	Eigen::Matrix<int,3,1> GetNCells();
93
+
94
+	/// Return a vector of the origin
95
+	Vector3r GetOrigin();
96
+
97
+	/// Return the centroid matrix
98
+	Vector3Xr GetCentroids();
99
+
100
+	/// Return a particular centroid
101
+	Vector3r GetCentroid(const int& vecpos);
102
+
103
+
104
+
105
+
106
+protected:
107
+
108
+	// ====================  LIFECYCLE     =======================
109
+
110
+	/// Default protected constructor. 
111
+	StructuredGrid ();                    
112
+
113
+	/// Default protected constructor. 
114
+	~StructuredGrid ();                    
115
+
116
+	// ====================  DATA MEMBERS  =========================
117
+
118
+	/// Number of cells in x,y,z
119
+	Eigen::Matrix<int,3,1> Nx3;
120
+
121
+	/// Cell spacings in each dimension (3 x n)
122
+	Vector3Xr Dx3;
123
+
124
+	/// Origin (3 x 1)
125
+	Vector3r Origin;
126
+
127
+	/// Vectorized edge coordinates
128
+	MatrixXr VecEdgCoord;
129
+
130
+	/// Vectorized centroid coordinates (3 x n)
131
+	Vector3Xr VecCentroid;
132
+
133
+	/// Log object
134
+	Logging* LogObj;
135
+
136
+	/// Vector of cell spacings in each dimension (3 x n)
137
+	std::vector<VectorXr> Dx3v;
138
+
139
+private:
140
+
141
+
142
+}; // -----  end of class  StructuredGrid  -----
143
+} // end of namespace formalhaut
144
+#endif

+ 47
- 0
src/exceptions.cpp View File

1
+// ===========================================================================
2
+// 
3
+//       Filename:  exceptions.cpp
4
+// 
5
+//    Description:  Holds the exceptions for the project
6
+// 
7
+//        Version:  0.0
8
+//        Created:  03/19/2013 12:18:41 PM
9
+//       Revision:  none
10
+//       Compiler:  Tested with g++ 
11
+// 
12
+//         Author:  M. Andy Kass (MAK)
13
+//
14
+//   Organisation:  US Geological Survey
15
+//                  
16
+//
17
+//          Email:  mkass@usgs.gov
18
+// 
19
+//  This program is free software: you can redistribute it and/or modify
20
+//  it under the terms of the GNU General Public License as published by
21
+//  the Free Software Foundation, either version 3 of the License, or
22
+//  (at your option) any later version.
23
+//
24
+//  This program is distributed in the hope that it will be useful,
25
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
26
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27
+//  GNU General Public License for more details.
28
+//
29
+//  You should have received a copy of the GNU General Public License
30
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
31
+//
32
+// ===========================================================================
33
+
34
+#include "exceptions.h"
35
+
36
+namespace formalhaut {
37
+
38
+	GenericFileIOError::GenericFileIOError(const std::string &filename) :
39
+			runtime_error("FILE I/O ERROR"){
40
+				std::cout << std::endl;
41
+				std::cout << "File I/O error while accessing " << 
42
+					filename << "." << std::endl;
43
+				std::cout << "Terminating." << std::endl;
44
+			}
45
+
46
+} // End of namespace formalhaut
47
+

+ 55
- 0
src/grid.cpp View File

1
+// ===========================================================================
2
+// 
3
+//       Filename:  grid.cpp
4
+// 
5
+//    Description:  
6
+// 
7
+//        Version:  0.0
8
+//        Created:  10/25/2009 07:42:08 AM
9
+//       Revision:  none
10
+//       Compiler:  g++ (c++)
11
+// 
12
+//         Author:  M. Andy Kass (MAK)
13
+//                  Trevor Irons (ti)
14
+//
15
+//   Organisation:  United States Geological Survey (USGS)
16
+//
17
+//          Email:  mkass@usgs.gov
18
+// 
19
+//  This program is free software: you can redistribute it and/or modify
20
+//  it under the terms of the GNU General Public License as published by
21
+//  the Free Software Foundation, either version 3 of the License, or
22
+//  (at your option) any later version.
23
+//
24
+//  This program is distributed in the hope that it will be useful,
25
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
26
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27
+//  GNU General Public License for more details.
28
+//
29
+//  You should have received a copy of the GNU General Public License
30
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
31
+//
32
+// ===========================================================================
33
+
34
+#include "grid.h"
35
+
36
+namespace formalhaut {
37
+
38
+	std::ostream &operator<<(std::ostream &stream, 
39
+				const Grid &ob) {
40
+		stream << "Class: Grid\n";
41
+		//stream << *(ReferenceCountedObject*)(&ob);
42
+		return stream;
43
+	}
44
+
45
+	Grid::Grid  ( ) {
46
+	}
47
+	
48
+	Grid::~Grid  ( ) {
49
+	}
50
+
51
+	gridtype Grid::GetGridType() {
52
+		return this->TypeOfGrid;
53
+	}
54
+
55
+}

+ 128
- 0
src/logging.cpp View File

1
+// ===========================================================================
2
+// 
3
+//       Filename:  logging.cpp
4
+// 
5
+//    Description:  
6
+// 
7
+//        Version:  0.0
8
+//        Created:  09/12/2013 03:26:27 PM
9
+//       Revision:  none
10
+//       Compiler:  Tested with g++ 
11
+// 
12
+//         Author:  M. Andy Kass (MAK)
13
+//
14
+//   Organisation:  US Geological Survey
15
+//                  
16
+//
17
+//          Email:  mkass@usgs.gov
18
+// 
19
+//  This program is free software: you can redistribute it and/or modify
20
+//  it under the terms of the GNU General Public License as published by
21
+//  the Free Software Foundation, either version 3 of the License, or
22
+//  (at your option) any later version.
23
+//
24
+//  This program is distributed in the hope that it will be useful,
25
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
26
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27
+//  GNU General Public License for more details.
28
+//
29
+//  You should have received a copy of the GNU General Public License
30
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
31
+//
32
+// ===========================================================================
33
+
34
+#include "logging.h"
35
+
36
+namespace formalhaut {
37
+
38
+	std::ostream &operator<<(std::ostream &stream,
39
+		const Logging &ob) {
40
+		stream << "Class: Logging\n";
41
+		stream << "Output file: " << ob.LogFile << "\n";
42
+		stream << "\n";
43
+
44
+
45
+		return stream;
46
+	}
47
+
48
+	Logging* Logging::New() {
49
+		return new Logging;
50
+	}
51
+
52
+	void Logging::Delete() {
53
+		if (this->outfile.is_open()) {
54
+			this->outfile.close();
55
+		}
56
+		delete this;
57
+	}
58
+
59
+	Logging::Logging() {
60
+
61
+	}
62
+
63
+	Logging::~Logging() {
64
+
65
+	}
66
+
67
+	void Logging::CloseLogFile() {
68
+
69
+
70
+		if (this->outfile.is_open()) {
71
+			this->outfile.close();
72
+		}
73
+
74
+	}
75
+
76
+	void Logging::WriteLogLine(const std::string line) {
77
+		
78
+		this->outfile << line << std::endl;
79
+
80
+	}
81
+
82
+	void Logging::WriteHeader(const std::string& name, const 
83
+		std::string& ver) {
84
+		
85
+		this->outfile << name << std::endl;
86
+		this->outfile << "version " << ver << std::endl << std::endl;
87
+
88
+		this->outfile << "Run began at " << this->GetDateTime() 
89
+			<< std::endl << std::endl;
90
+	}
91
+
92
+	std::string Logging::GetDateTime() {
93
+		time_t now = time(0);
94
+		struct tm tstruct;
95
+		char buf[80];
96
+		tstruct = *localtime(&now);
97
+		strftime(buf,sizeof(buf),"%Y-%m-%d.%X",&tstruct);
98
+		return buf;
99
+	
100
+	}
101
+
102
+	void Logging::WriteCloser() {
103
+
104
+		this->outfile << "Successfully ended at " << this->GetDateTime()
105
+			<< std::endl;
106
+		this->outfile.close();
107
+	}
108
+
109
+	void Logging::WriteDashedLine() {
110
+		std::string dashes;
111
+		dashes = "------------------------------------";
112
+		dashes = dashes + dashes;
113
+		this->outfile << dashes << std::endl;
114
+	}
115
+
116
+	void Logging::SetLogFile(const std::string& fname) {
117
+
118
+		this->LogFile = fname;
119
+		this->outfile.open(this->LogFile.c_str());
120
+	}
121
+
122
+	std::string Logging::GetLogFile() {
123
+		return this->LogFile;
124
+	}
125
+
126
+
127
+
128
+} // end of namespace formalhaut

+ 433
- 0
src/rectilinearstructuredgrid.cpp View File

1
+// ===========================================================================
2
+// 
3
+//       Filename:  rectilinearstructuredgrid.cpp
4
+// 
5
+//    Description:  
6
+// 
7
+//        Version:  0.0
8
+//        Created:  09/16/2013 01:55:17 PM
9
+//       Revision:  none
10
+//       Compiler:  Tested with g++ 
11
+// 
12
+//         Author:  M. Andy Kass (MAK)
13
+//
14
+//   Organisation:  US Geological Survey
15
+//                  
16
+//
17
+//          Email:  mkass@usgs.gov
18
+// 
19
+//  This program is free software: you can redistribute it and/or modify
20
+//  it under the terms of the GNU General Public License as published by
21
+//  the Free Software Foundation, either version 3 of the License, or
22
+//  (at your option) any later version.
23
+//
24
+//  This program is distributed in the hope that it will be useful,
25
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
26
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27
+//  GNU General Public License for more details.
28
+//
29
+//  You should have received a copy of the GNU General Public License
30
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
31
+//
32
+// ===========================================================================
33
+
34
+#include "rectilinearstructuredgrid.h"
35
+
36
+namespace formalhaut {
37
+
38
+	std::ostream &operator<<(std::ostream &stream,
39
+		const RectilinearStructuredGrid &ob) {
40
+
41
+		return stream;
42
+	}
43
+
44
+	RectilinearStructuredGrid::RectilinearStructuredGrid() {
45
+		this->TypeOfGrid = rectilinear;
46
+	}
47
+
48
+	RectilinearStructuredGrid::~RectilinearStructuredGrid() {
49
+
50
+	}
51
+
52
+	RectilinearStructuredGrid* RectilinearStructuredGrid::New() {
53
+		return new RectilinearStructuredGrid;
54
+	}
55
+
56
+	void RectilinearStructuredGrid::Delete() {
57
+		delete this;
58
+	}
59
+
60
+	void RectilinearStructuredGrid::ComputeCentroids() {
61
+
62
+		int ntot;
63
+		ntot = this->Nx3(0)*this->Nx3(1)*this->Nx3(2);
64
+		this->VecCentroid.resize(Eigen::NoChange,ntot);
65
+		this->CW.resize(Eigen::NoChange,ntot);
66
+		int mm = 0;
67
+
68
+		// First z changes, then x then y
69
+		// first depth, then easting, then northing
70
+		//Original
71
+		
72
+		///*
73
+		//std::cout << this->Dx3v[0].segment(0,1).sum() << std::endl;
74
+		//std::cout << this->Dx3v[0](0) << "  " << this->Dx3v[0](1) <<
75
+		//	std::endl;
76
+		//std::cout << this->Dx3v[0].segment(0,1).sum() - 
77
+		//	this->Dx3v[0](1)/2 + this->Dx3v[0](0) + 
78
+		//	this->Origin(0) << std::endl;
79
+		for (int jj=0;jj<this->Nx3(1);jj++) {
80
+			for (int ii=this->Nx3(0)-1;ii>-1;ii--) {
81
+			//for (int ii=0;ii<this->Nx3(0);ii++) {
82
+				for (int kk=0;kk<this->Nx3(2);kk++) {
83
+		
84
+		//*/
85
+		/*
86
+		for (int ii=0;ii<this->Nx3(0);ii++) {
87
+			for (int jj=0;jj<this->Nx3(1);jj++) {
88
+				for (int kk=0;kk<this->Nx3(2);kk++) {
89
+		*/
90
+					this->VecCentroid(0,mm) = this->Origin(0) +
91
+						this->Dx3v[0].segment(0,ii).sum() - 
92
+						this->Dx3v[0](ii)/2 + this->Dx3v[0](ii);
93
+					this->VecCentroid(1,mm) = this->Origin(1) +
94
+						this->Dx3v[1].segment(0,jj).sum() -
95
+						this->Dx3v[1](jj)/2 + this->Dx3v[1](jj);
96
+					this->VecCentroid(2,mm) = this->Origin(2) +
97
+						this->Dx3v[2].segment(0,kk).sum() -
98
+						this->Dx3v[2](kk)/2 + this->Dx3v[2](kk);
99
+					this->CW(0,mm) = this->Dx3v[0](ii);
100
+					this->CW(1,mm) = this->Dx3v[1](jj);
101
+					this->CW(2,mm) = this->Dx3v[2](kk);
102
+					mm++;
103
+				}
104
+			}
105
+		}
106
+
107
+	}
108
+
109
+	void RectilinearStructuredGrid::SetCellDims(const
110
+		std::vector<VectorXr>& dx3v) {
111
+
112
+		this->Dx3v[0] = dx3v[0];
113
+		this->Dx3v[1] = dx3v[1];
114
+		this->Dx3v[2] = dx3v[2];
115
+
116
+	}
117
+
118
+	std::vector<vtkActor*> RectilinearStructuredGrid::ReturnMeshActor() {
119
+
120
+		// Coords are edges of cells, not centroids
121
+		vtkFloatArray *xCoords = vtkFloatArray::New();
122
+		vtkFloatArray *yCoords = vtkFloatArray::New();
123
+		vtkFloatArray *zCoords = vtkFloatArray::New();
124
+		Real val;
125
+
126
+		val = this->Origin(0);
127
+		xCoords->InsertNextValue(val);
128
+		for (int ii=1;ii<=this->Nx3(0)+1;ii++) {
129
+			val = this->Origin(0) + this->Dx3v[0].segment(0,ii-1).sum();
130
+			xCoords->InsertNextValue(val);
131
+		}
132
+		val = this->Origin(1);
133
+		yCoords->InsertNextValue(val);
134
+		for (int ii=1;ii<=this->Nx3(1)+1;ii++) {
135
+			val = this->Origin(1) + this->Dx3v[1].segment(0,ii-1).sum();
136
+			yCoords->InsertNextValue(val);
137
+		}
138
+		val = -1 * this->Origin(2);
139
+		zCoords->InsertNextValue(val);
140
+		for (int ii=1;ii<=this->Nx3(2)+1;ii++) {
141
+			val = this->Origin(2) + this->Dx3v[2].segment(0,ii-1).sum();
142
+			val = -1*val;
143
+			zCoords->InsertNextValue(val);
144
+		}
145
+
146
+		vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
147
+		rgrid->SetDimensions(this->Nx3(1)+2,this->Nx3(0)+2,this->Nx3(2)+2);
148
+		rgrid->SetXCoordinates(yCoords);
149
+		rgrid->SetYCoordinates(xCoords);
150
+		rgrid->SetZCoordinates(zCoords);
151
+
152
+		vtkRectilinearGridGeometryFilter *plane = 
153
+			vtkRectilinearGridGeometryFilter::New();
154
+		plane->SetInputData(rgrid);
155
+		//plane->SetExtent(0,this->Nx3(0),0,this->Nx3(1),0,this->Nx3(2));
156
+		//Just for testing
157
+		plane->SetExtent(0,0,0,this->Nx3(0)+1,0,this->Nx3(2)+1);
158
+
159
+		vtkPolyDataMapper *rgridmapper = vtkPolyDataMapper::New();
160
+		rgridmapper->SetInputConnection(plane->GetOutputPort());
161
+
162
+		vtkActor *wireActor = vtkActor::New();
163
+		wireActor->SetMapper(rgridmapper);
164
+		wireActor->GetProperty()->SetRepresentationToWireframe();
165
+		wireActor->GetProperty()->SetColor(0,1,0);
166
+
167
+		//Make an actor for each face
168
+		//Second face
169
+		vtkRectilinearGridGeometryFilter *plane1 =
170
+			vtkRectilinearGridGeometryFilter::New();
171
+		plane1->SetInputData(rgrid);
172
+		plane1->SetExtent(this->Nx3(1)+1,this->Nx3(1)+1,0,this->Nx3(0)+1,
173
+			0,this->Nx3(2)+1);
174
+		vtkPolyDataMapper *rgridmapper1 = vtkPolyDataMapper::New();
175
+		rgridmapper1->SetInputConnection(plane1->GetOutputPort());
176
+		vtkActor *wireActor1 = vtkActor::New();
177
+		wireActor1->SetMapper(rgridmapper1);
178
+		wireActor1->GetProperty()->SetRepresentationToWireframe();
179
+		wireActor1->GetProperty()->SetColor(0,1,0);
180
+		//Third face
181
+		vtkRectilinearGridGeometryFilter *plane2 =
182
+			vtkRectilinearGridGeometryFilter::New();
183
+		plane2->SetInputData(rgrid);
184
+		plane2->SetExtent(0,this->Nx3(1)+1,0,0,
185
+			0,this->Nx3(2)+1);
186
+		vtkPolyDataMapper *rgridmapper2 = vtkPolyDataMapper::New();
187
+		rgridmapper2->SetInputConnection(plane2->GetOutputPort());
188
+		vtkActor *wireActor2 = vtkActor::New();
189
+		wireActor2->SetMapper(rgridmapper2);
190
+		wireActor2->GetProperty()->SetRepresentationToWireframe();
191
+		wireActor2->GetProperty()->SetColor(0,1,0);
192
+		//Fourth face
193
+		vtkRectilinearGridGeometryFilter *plane3 =
194
+			vtkRectilinearGridGeometryFilter::New();
195
+		plane3->SetInputData(rgrid);
196
+		plane3->SetExtent(0,this->Nx3(1)+1,this->Nx3(0)+1,this->Nx3(0)+1,
197
+			0,this->Nx3(2)+1);
198
+		vtkPolyDataMapper *rgridmapper3 = vtkPolyDataMapper::New();
199
+		rgridmapper3->SetInputConnection(plane3->GetOutputPort());
200
+		vtkActor *wireActor3 = vtkActor::New();
201
+		wireActor3->SetMapper(rgridmapper3);
202
+		wireActor3->GetProperty()->SetRepresentationToWireframe();
203
+		wireActor3->GetProperty()->SetColor(0,1,0);
204
+		//Fifth face
205
+		vtkRectilinearGridGeometryFilter *plane4 =
206
+			vtkRectilinearGridGeometryFilter::New();
207
+		plane4->SetInputData(rgrid);
208
+		plane4->SetExtent(0,this->Nx3(1)+1,0,this->Nx3(0)+1,
209
+			0,0);
210
+		vtkPolyDataMapper *rgridmapper4 = vtkPolyDataMapper::New();
211
+		rgridmapper4->SetInputConnection(plane4->GetOutputPort());
212
+		vtkActor *wireActor4 = vtkActor::New();
213
+		wireActor4->SetMapper(rgridmapper4);
214
+		wireActor4->GetProperty()->SetRepresentationToWireframe();
215
+		wireActor4->GetProperty()->SetColor(0,1,0);
216
+		//Sixth face
217
+		vtkRectilinearGridGeometryFilter *plane5 =
218
+			vtkRectilinearGridGeometryFilter::New();
219
+		plane5->SetInputData(rgrid);
220
+		plane5->SetExtent(0,this->Nx3(1)+1,0,this->Nx3(0)+1,
221
+			this->Nx3(2)+1,this->Nx3(2)+1);
222
+		vtkPolyDataMapper *rgridmapper5 = vtkPolyDataMapper::New();
223
+		rgridmapper5->SetInputConnection(plane5->GetOutputPort());
224
+		vtkActor *wireActor5 = vtkActor::New();
225
+		wireActor5->SetMapper(rgridmapper5);
226
+		wireActor5->GetProperty()->SetRepresentationToWireframe();
227
+		wireActor5->GetProperty()->SetColor(0,1,0);
228
+		
229
+		vtkSmartPointer<vtkOutlineFilter> outlineFilter = 
230
+			vtkSmartPointer<vtkOutlineFilter>::New();
231
+		outlineFilter->SetInputData(rgrid);
232
+		vtkSmartPointer<vtkDataSetMapper> outlinemapper = 
233
+			vtkSmartPointer<vtkDataSetMapper>::New();
234
+		outlinemapper->SetInputConnection(outlineFilter->GetOutputPort());
235
+
236
+		vtkActor *outlineActor = vtkActor::New();
237
+		outlineActor->SetMapper(outlinemapper);
238
+
239
+		std::vector<vtkActor*> actorreturn;
240
+		// Order of actors:
241
+		// 1. outlineActor
242
+		// 2-7. plane 0-5
243
+		
244
+		actorreturn.push_back(outlineActor);
245
+		actorreturn.push_back(wireActor);
246
+		actorreturn.push_back(wireActor1);
247
+		actorreturn.push_back(wireActor2);
248
+		actorreturn.push_back(wireActor3);
249
+		actorreturn.push_back(wireActor4);
250
+		actorreturn.push_back(wireActor5);
251
+
252
+		return actorreturn;
253
+
254
+
255
+	}
256
+
257
+	std::vector<VectorXr> RectilinearStructuredGrid::GetCellDims() {
258
+		return this->Dx3v;
259
+	}
260
+
261
+	VectorXr RectilinearStructuredGrid::GetCorners(const int& vecpos) {
262
+		VectorXr corners;
263
+		corners.resize(6);
264
+		Vector3r cntr;
265
+
266
+		cntr = this->VecCentroid.col(vecpos);
267
+
268
+		corners(0) = cntr(0) - this->CW(0,vecpos)/2;
269
+		corners(1) = cntr(0) + this->CW(0,vecpos)/2;
270
+		corners(2) = cntr(1) - this->CW(1,vecpos)/2;
271
+		corners(3) = cntr(1) + this->CW(1,vecpos)/2;
272
+		corners(4) = cntr(2) - this->CW(2,vecpos)/2;
273
+		corners(5) = cntr(2) + this->CW(2,vecpos)/2;
274
+		
275
+		return corners;
276
+	}
277
+
278
+	void RectilinearStructuredGrid::DisplayMesh() {
279
+
280
+		// Coords are edges of cells, not centroids
281
+		vtkFloatArray *xCoords = vtkFloatArray::New();
282
+		vtkFloatArray *yCoords = vtkFloatArray::New();
283
+		vtkFloatArray *zCoords = vtkFloatArray::New();
284
+		Real val;
285
+
286
+		val = this->Origin(0);
287
+		xCoords->InsertNextValue(val);
288
+		for (int ii=1;ii<=this->Nx3(0)+1;ii++) {
289
+			val = this->Origin(0) + this->Dx3v[0].segment(0,ii-1).sum();
290
+			xCoords->InsertNextValue(val);
291
+		}
292
+		val = this->Origin(1);
293
+		yCoords->InsertNextValue(val);
294
+		for (int ii=1;ii<=this->Nx3(1)+1;ii++) {
295
+			val = this->Origin(1) + this->Dx3v[1].segment(0,ii-1).sum();
296
+			yCoords->InsertNextValue(val);
297
+		}
298
+		val = -1 * this->Origin(2);
299
+		zCoords->InsertNextValue(val);
300
+		for (int ii=1;ii<=this->Nx3(2)+1;ii++) {
301
+			val = this->Origin(2) + this->Dx3v[2].segment(0,ii-1).sum();
302
+			val = -1*val;
303
+			zCoords->InsertNextValue(val);
304
+		}
305
+
306
+		vtkRectilinearGrid *rgrid = vtkRectilinearGrid::New();
307
+		rgrid->SetDimensions(this->Nx3(1)+2,this->Nx3(0)+2,this->Nx3(2)+2);
308
+		rgrid->SetXCoordinates(yCoords);
309
+		rgrid->SetYCoordinates(xCoords);
310
+		rgrid->SetZCoordinates(zCoords);
311
+
312
+		vtkRectilinearGridGeometryFilter *plane = 
313
+			vtkRectilinearGridGeometryFilter::New();
314
+		plane->SetInputData(rgrid);
315
+		//plane->SetExtent(0,this->Nx3(0),0,this->Nx3(1),0,this->Nx3(2));
316
+		//Just for testing
317
+		plane->SetExtent(0,0,0,this->Nx3(0)+1,0,this->Nx3(2)+1);
318
+
319
+		vtkPolyDataMapper *rgridmapper = vtkPolyDataMapper::New();
320
+		rgridmapper->SetInputConnection(plane->GetOutputPort());
321
+
322
+		vtkActor *wireActor = vtkActor::New();
323
+		wireActor->SetMapper(rgridmapper);
324
+		wireActor->GetProperty()->SetRepresentationToWireframe();
325
+		wireActor->GetProperty()->SetColor(0,1,0);
326
+
327
+		//Make an actor for each face
328
+		//Second face
329
+		vtkRectilinearGridGeometryFilter *plane1 =
330
+			vtkRectilinearGridGeometryFilter::New();
331
+		plane1->SetInputData(rgrid);
332
+		plane1->SetExtent(this->Nx3(1)+1,this->Nx3(1)+1,0,this->Nx3(0)+1,
333
+			0,this->Nx3(2)+1);
334
+		vtkPolyDataMapper *rgridmapper1 = vtkPolyDataMapper::New();
335
+		rgridmapper1->SetInputConnection(plane1->GetOutputPort());
336
+		vtkActor *wireActor1 = vtkActor::New();
337
+		wireActor1->SetMapper(rgridmapper1);
338
+		wireActor1->GetProperty()->SetRepresentationToWireframe();
339
+		wireActor1->GetProperty()->SetColor(0,1,0);
340
+		//Third face
341
+		vtkRectilinearGridGeometryFilter *plane2 =
342
+			vtkRectilinearGridGeometryFilter::New();
343
+		plane2->SetInputData(rgrid);
344
+		plane2->SetExtent(0,this->Nx3(1)+1,0,0,
345
+			0,this->Nx3(2)+1);
346
+		vtkPolyDataMapper *rgridmapper2 = vtkPolyDataMapper::New();
347
+		rgridmapper2->SetInputConnection(plane2->GetOutputPort());
348
+		vtkActor *wireActor2 = vtkActor::New();
349
+		wireActor2->SetMapper(rgridmapper2);
350
+		wireActor2->GetProperty()->SetRepresentationToWireframe();
351
+		wireActor2->GetProperty()->SetColor(0,1,0);
352
+		//Fourth face
353
+		vtkRectilinearGridGeometryFilter *plane3 =
354
+			vtkRectilinearGridGeometryFilter::New();
355
+		plane3->SetInputData(rgrid);
356
+		plane3->SetExtent(0,this->Nx3(1)+1,this->Nx3(0)+1,this->Nx3(0)+1,
357
+			0,this->Nx3(2)+1);
358
+		vtkPolyDataMapper *rgridmapper3 = vtkPolyDataMapper::New();
359
+		rgridmapper3->SetInputConnection(plane3->GetOutputPort());
360
+		vtkActor *wireActor3 = vtkActor::New();
361
+		wireActor3->SetMapper(rgridmapper3);
362
+		wireActor3->GetProperty()->SetRepresentationToWireframe();
363
+		wireActor3->GetProperty()->SetColor(0,1,0);
364
+		//Fifth face
365
+		vtkRectilinearGridGeometryFilter *plane4 =
366
+			vtkRectilinearGridGeometryFilter::New();
367
+		plane4->SetInputData(rgrid);
368
+		plane4->SetExtent(0,this->Nx3(1)+1,0,this->Nx3(0)+1,
369
+			0,0);
370
+		vtkPolyDataMapper *rgridmapper4 = vtkPolyDataMapper::New();
371
+		rgridmapper4->SetInputConnection(plane4->GetOutputPort());
372
+		vtkActor *wireActor4 = vtkActor::New();
373
+		wireActor4->SetMapper(rgridmapper4);
374
+		wireActor4->GetProperty()->SetRepresentationToWireframe();
375
+		wireActor4->GetProperty()->SetColor(0,1,0);
376
+		//Sixth face
377
+		vtkRectilinearGridGeometryFilter *plane5 =
378
+			vtkRectilinearGridGeometryFilter::New();
379
+		plane5->SetInputData(rgrid);
380
+		plane5->SetExtent(0,this->Nx3(1)+1,0,this->Nx3(0)+1,
381
+			this->Nx3(2)+1,this->Nx3(2)+1);
382
+		vtkPolyDataMapper *rgridmapper5 = vtkPolyDataMapper::New();
383
+		rgridmapper5->SetInputConnection(plane5->GetOutputPort());
384
+		vtkActor *wireActor5 = vtkActor::New();
385
+		wireActor5->SetMapper(rgridmapper5);
386
+		wireActor5->GetProperty()->SetRepresentationToWireframe();
387
+		wireActor5->GetProperty()->SetColor(0,1,0);
388
+		
389
+		vtkSmartPointer<vtkOutlineFilter> outlineFilter = 
390
+			vtkSmartPointer<vtkOutlineFilter>::New();
391
+		outlineFilter->SetInputData(rgrid);
392
+		vtkSmartPointer<vtkDataSetMapper> outlinemapper = 
393
+			vtkSmartPointer<vtkDataSetMapper>::New();
394
+		outlinemapper->SetInputConnection(outlineFilter->GetOutputPort());
395
+
396
+		vtkActor *outlineActor = vtkActor::New();
397
+		outlineActor->SetMapper(outlinemapper);
398
+
399
+		vtkRenderer *renderer = vtkRenderer::New();
400
+		vtkRenderWindow *renWin = vtkRenderWindow::New();
401
+		renWin->AddRenderer(renderer);
402
+		vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
403
+		iren->SetRenderWindow(renWin);
404
+
405
+		renderer->AddActor(wireActor);
406
+		//renderer->AddActor(wireActor1);
407
+		renderer->AddActor(wireActor2);
408
+		//renderer->AddActor(wireActor3);
409
+		//renderer->AddActor(wireActor4);
410
+		renderer->AddActor(wireActor5);
411
+		renderer->AddActor(outlineActor);
412
+		renderer->SetBackground(0,0,0);
413
+		renderer->ResetCamera();
414
+		renderer->GetActiveCamera()->Elevation(40.0);
415
+		renderer->GetActiveCamera()->Azimuth(15);
416
+		renderer->GetActiveCamera()->Roll(200);
417
+		renderer->GetActiveCamera()->Zoom(1.0);
418
+
419
+		renWin->SetSize(900,900);
420
+
421
+		renWin->Render();
422
+		iren->Start();
423
+
424
+
425
+
426
+
427
+
428
+	}
429
+
430
+
431
+
432
+
433
+}

+ 470
- 0
src/structuredgrid.cpp View File

1
+// ===========================================================================
2
+// 
3
+//       Filename:  structuredgrid.cpp
4
+// 
5
+//    Description:  
6
+// 
7
+//        Version:  0.0
8
+//        Created:  09/12/2013 08:36:55 AM
9
+//       Revision:  none
10
+//       Compiler:  Tested with g++ 
11
+// 
12
+//         Author:  M. Andy Kass (MAK)
13
+//
14
+//   Organisation:  US Geological Survey
15
+//                  
16
+//
17
+//          Email:  mkass@usgs.gov
18
+// 
19
+//  This program is free software: you can redistribute it and/or modify
20
+//  it under the terms of the GNU General Public License as published by
21
+//  the Free Software Foundation, either version 3 of the License, or
22
+//  (at your option) any later version.
23
+//
24
+//  This program is distributed in the hope that it will be useful,
25
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
26
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27
+//  GNU General Public License for more details.
28
+//
29
+//  You should have received a copy of the GNU General Public License
30
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
31
+//
32
+// ===========================================================================
33
+
34
+#include "structuredgrid.h"
35
+
36
+namespace formalhaut {
37
+
38
+	std::ostream &operator<<(std::ostream &stream,
39
+		const StructuredGrid &ob) {
40
+		stream << "Class: StructuredGrid\n";
41
+
42
+		return stream;
43
+	}
44
+
45
+	StructuredGrid::StructuredGrid() : LogObj(NULL) {
46
+
47
+	}
48
+
49
+	StructuredGrid::~StructuredGrid() {
50
+
51
+	}
52
+
53
+	void StructuredGrid::SetNCells(const Eigen::Matrix<int,3,1>& ncells) {
54
+		this->Nx3 = ncells;
55
+	}
56
+
57
+	void StructuredGrid::SetOrigin(const Vector3r& origin) {
58
+		this->Origin = origin;
59
+	}
60
+
61
+	void StructuredGrid::ReadMeshtoolsGrid(const std::string& fname) {
62
+
63
+		Eigen::Matrix<int,3,1> tn;
64
+		Vector3r to;
65
+		VectorXr dx;
66
+		VectorXr dy;
67
+		VectorXr dz;
68
+
69
+		int adds = 0;
70
+		int multiplier;
71
+		int bufsize;
72
+		std::string temp;
73
+		std::string temp2;
74
+		std::string tast;
75
+		char* commentBuffer = new char[200];
76
+
77
+		bufsize = 200;
78
+		std::fstream infile(fname.c_str(),std::ios::in);
79
+		if (infile.fail()) {
80
+			throw GenericFileIOError(fname);
81
+		}
82
+		infile.getline(commentBuffer,bufsize);
83
+		temp = commentBuffer;
84
+		std::istringstream in(temp);
85
+		for (int ii=0;ii<3;ii++) {
86
+			in >> tn(ii);
87
+		}
88
+		in.clear();
89
+
90
+		infile.getline(commentBuffer,bufsize);
91
+		temp = " ";
92
+		temp = commentBuffer;
93
+		std::istringstream in2(temp);
94
+		for (int ii=0;ii<=2;ii++) {
95
+			in2 >> to(ii);
96
+		}
97
+		in2.clear();
98
+
99
+		dx.resize(tn(0));
100
+		dy.resize(tn(1));
101
+		dz.resize(tn(2));
102
+
103
+		infile.getline(commentBuffer,bufsize);
104
+		temp = commentBuffer;
105
+		in.str(temp.c_str());
106
+		int ii = 0;
107
+		Real t2;
108
+		while (ii <=tn(0)-1) {
109
+			in >> t2;
110
+			tast = in.peek();
111
+			if (tast=="*") {
112
+				temp2 = in.get();
113
+				multiplier = t2;
114
+				in >> t2;
115
+				for (int jj=0;jj<=multiplier-2;jj++) {
116
+					dx(ii+jj) = t2;
117
+					adds = jj;
118
+				}
119
+				ii = ii+adds+1;
120
+			}
121
+			dx(ii) = t2;
122
+			ii++;
123
+		}
124
+
125
+		infile.getline(commentBuffer,bufsize);
126
+		temp = commentBuffer;
127
+		in.str(temp.c_str());
128
+		ii = 0;
129
+		adds = 0;
130
+		while (ii <=tn(1)-1) {
131
+			in >> t2;
132
+			tast = in.peek();
133
+			if (tast=="*") {
134
+				temp2 = in.get();
135
+				multiplier = t2;
136
+				in >> t2;
137
+				for (int jj=0;jj<=multiplier-2;jj++) {
138
+					dy(ii+jj) = t2;
139
+					adds = jj;
140
+				}
141
+				ii = ii+adds+1;
142
+			}
143
+			dy(ii) = t2;
144
+			ii++;
145
+		}
146
+
147
+		infile.getline(commentBuffer,bufsize);
148
+		temp = commentBuffer;
149
+		in.str(temp.c_str());
150
+		ii = 0;
151
+		adds = 0;
152
+		while (ii <=tn(2)-1) {
153
+			in >> t2;
154
+			tast = in.peek();
155
+			if (tast=="*") {
156
+				temp2 = in.get();
157
+				multiplier = t2;
158
+				in >> t2;
159
+				for (int jj=0;jj<=multiplier-2;jj++) {
160
+					dz(ii+jj) = t2;
161
+					adds = jj;
162
+				}
163
+				ii = ii+adds+1;
164
+			}
165
+			dz(ii) = t2;
166
+			ii++;
167
+		}
168
+
169
+		
170
+
171
+/*		
172
+		infile.getline(commentBuffer,bufsize);
173
+		temp = commentBuffer;
174
+		in.str(temp.c_str());
175
+		int ii = 0;
176
+		while (ii <= tn(0)-1) {
177
+			in >> dx(ii);
178
+			std::cout << dx(ii) << std::endl;
179
+			ii++;
180
+			tast = in.peek();
181
+			if (tast=="*") {
182
+				temp2 = in.get();
183
+				in >> multiplier;
184
+				std::cout << multiplier << std::endl;
185
+				for (int jj=0;jj<=multiplier-2;jj++) {
186
+					dx(ii+jj) = dx(ii-1);
187
+					adds = jj;
188
+				}
189
+				ii = ii+adds+1;
190
+			}
191
+		}
192
+		in.clear();
193
+
194
+		infile.getline(commentBuffer,bufsize);
195
+		temp = commentBuffer;
196
+		in.str(temp.c_str());
197
+		ii = 0;
198
+		adds = 0;
199
+		while (ii <= tn(1)-1) {
200
+			in >> dy(ii);
201
+			ii++;
202
+			tast = in.peek();
203
+			if (tast == "*") {
204
+				temp2 = in.get();
205
+				in >> multiplier;
206
+				for (int jj=0;jj<=multiplier-2;jj++) {
207
+					dy(ii+jj) = dy(ii-1);
208
+					adds = jj;
209
+				}
210
+				ii = ii+adds+1;
211
+			}
212
+		}
213
+		in.clear();
214
+
215
+		infile.getline(commentBuffer,bufsize);
216
+		temp = commentBuffer;
217
+		in.str(temp.c_str());
218
+		ii = 0;
219
+		adds = 0;
220
+		while (ii <= tn(2)-1) {
221
+			in >> dz(ii);
222
+			ii++;
223
+			tast = in.peek();
224
+			if (tast=="*") {
225
+				temp2=in.get();
226
+				in >> multiplier;
227
+				for (int jj=0;jj<=multiplier-2;jj++) {
228
+					dz(ii+jj) = dz(ii-1);
229
+					adds = jj;
230
+				}
231
+				ii = ii+adds+1;
232
+			}
233
+		}
234
+		in.clear();
235
+*/
236
+		infile.close();
237
+		delete [] commentBuffer;
238
+
239
+		this->Origin = to;
240
+		this->Nx3 = tn;
241
+
242
+		if (this->TypeOfGrid == regular) {
243
+			// Check to make sure mesh is regular
244
+			for (int jj=0;jj<this->Nx3(0);jj++) {
245
+				if (dx(0) != dx(jj)) {
246
+					std::cout << "Error!  Non-regular mesh read into " <<
247
+						"regular mesh object.  Exiting." << std::endl;
248
+					exit(EXIT_FAILURE);
249
+				}
250
+			}
251
+			for (int jj=0;jj<this->Nx3(1);jj++) {
252
+				if (dy(0) != dy(jj)) {
253
+					std::cout << "Error!  Non-regular mesh read into " <<
254
+						"regular mesh object.  Exiting." << std::endl;
255
+					exit(EXIT_FAILURE);
256
+				}
257
+			}
258
+			for (int jj=0;jj<this->Nx3(2);jj++) {
259
+				if (dz(0) != dz(jj)) {
260
+					std::cout << "Error!  Non-regular mesh read into " <<
261
+						"regular mesh object.  Exiting." << std::endl;
262
+					exit(EXIT_FAILURE);
263
+				}
264
+			}
265
+
266
+			// Resize Dx3 and populate
267
+			this->Dx3.resize(Eigen::NoChange,1);
268
+			this->Dx3(0) = dx(0);
269
+			this->Dx3(1) = dy(0);
270
+			this->Dx3(2) = dz(0);
271
+
272
+			// Write to log file
273
+			std::string logw;
274
+			std::stringstream ss;
275
+			std::vector<std::string> logs;
276
+			this->LogObj->WriteLogLine("Regular Mesh parameters:");
277
+			logw = "nx\t\tny\t\tnz\t\tdx\t\tdy\t\tdz\t\tox\t\toy\t\toz";
278
+			this->LogObj->WriteLogLine(logw);
279
+			this->LogObj->WriteDashedLine();
280
+			ss << this->Nx3(0);
281
+			logs.push_back(ss.str());
282
+			ss.str("");
283
+			ss << this->Nx3(1);
284
+			logs.push_back(ss.str());
285
+			ss.str("");
286
+			ss << this->Nx3(2);
287
+			logs.push_back(ss.str());
288
+			ss.str("");
289
+			ss << this->Dx3(0);
290
+			logs.push_back(ss.str());
291
+			ss.str("");
292
+			ss << this->Dx3(1);
293
+			logs.push_back(ss.str());
294
+			ss.str("");
295
+			ss << this->Dx3(2);
296
+			logs.push_back(ss.str());
297
+			ss.str("");
298
+			ss << this->Origin(0);
299
+			logs.push_back(ss.str());
300
+			ss.str("");
301
+			ss << this->Origin(1);
302
+			logs.push_back(ss.str());
303
+			ss.str("");
304
+			ss << this->Origin(2);
305
+			logs.push_back(ss.str());
306
+			ss.str("");
307
+			logw = "";
308
+			for (int jj=0;jj<9;jj++) {
309
+				logw = logw + logs[jj] + "\t\t";
310
+			}
311
+			this->LogObj->WriteLogLine(logw);
312
+			this->LogObj->WriteLogLine("\n");
313
+			
314
+
315
+		}
316
+		else if (this->TypeOfGrid == rectilinear) {
317
+			this->Dx3v.push_back(dx);
318
+			this->Dx3v.push_back(dy);
319
+			this->Dx3v.push_back(dz);
320
+
321
+			// Write to log file
322
+			std::string logw;
323
+			std::stringstream ss;
324
+			std::vector<std::string> logs;
325
+			this->LogObj->WriteLogLine("Rectilinear Mesh Parameters:");
326
+			logw = "nx\t\tny\t\tnz\t\tox\t\toy\t\toz";
327
+			this->LogObj->WriteLogLine(logw);
328
+			this->LogObj->WriteDashedLine();
329
+			ss << this->Nx3(0);
330
+			logs.push_back(ss.str());
331
+			ss.str("");
332
+			ss << this->Nx3(1);
333
+			logs.push_back(ss.str());
334
+			ss.str("");
335
+			ss << this->Nx3(2);
336
+			logs.push_back(ss.str());
337
+			ss.str("");
338
+			ss << this->Origin(0);
339
+			logs.push_back(ss.str());
340
+			ss.str("");
341
+			ss << this->Origin(1);
342
+			logs.push_back(ss.str());
343
+			ss.str("");
344
+			ss << this->Origin(2);
345
+			logs.push_back(ss.str());
346
+			ss.str("");
347
+			logw = "";
348
+			for (int jj=0;jj<6;jj++) {
349
+				logw = logw + logs[jj] + "\t\t";
350
+			}
351
+			this->LogObj->WriteLogLine(logw);
352
+			this->LogObj->WriteLogLine("\n");
353
+
354
+
355
+		}
356
+		else if (this->TypeOfGrid == octree) {
357
+			std::cout << "Error!  Octree mesh object cannot interface " <<
358
+				"with Meshtools format.  Exiting." << std::endl;
359
+			exit(EXIT_FAILURE);
360
+		}
361
+		else {
362
+			std::cout << "Error!  Mesh object not declared properly. Exiting."
363
+				<< std::endl << "(will try to be more unhelpful in the "
364
+				<< "future)" << std::endl;
365
+			exit(EXIT_FAILURE);
366
+		}
367
+
368
+
369
+	}
370
+
371
+	void StructuredGrid::AttachLogObject(Logging* logobj) {
372
+		this->LogObj = logobj;
373
+	}
374
+
375
+	int StructuredGrid::GetVectorPosition(const Eigen::Matrix<int,3,1>& 
376
+		cellpos) {
377
+
378
+		int vecpos;
379
+
380
+		// Check and make sure the points are valid
381
+		if (cellpos(0) >= this->Nx3(0)) {
382
+			std::cout << "Invalid queried vector position in x dim\n";
383
+			exit(EXIT_FAILURE);
384
+		}
385
+		if (cellpos(1) >= this->Nx3(1)) {
386
+			std::cout << "Invalid queried vector position in y dim\n";
387
+			exit(EXIT_FAILURE);
388
+		}
389
+		if (cellpos(2) >= this->Nx3(2)) {
390
+			std::cout << "Invalid queried vector position in z dim\n";
391
+			exit(EXIT_FAILURE);
392
+		}
393
+
394
+		// z changes first, then x, then y
395
+		vecpos = (cellpos(1)*this->Nx3(0)*this->Nx3(2)) +
396
+			(cellpos(0)*this->Nx3(2)) + cellpos(2);
397
+
398
+		return vecpos;
399
+	}
400
+
401
+	int StructuredGrid::QueryCellPosition(const Vector3r& cellpos) {
402
+		int vecpos;
403
+		Vector3r cellposmo;
404
+
405
+		//Check to make sure the point is in the model region
406
+		if (cellpos(0) < this->Origin(0) || cellpos(0) >= 
407
+			this->Nx3(0)*this->Dx3(0) + this->Origin(0)) {
408
+			std::cout << "Queried point exists outside model region.\n";
409
+			std::cout << "Given location: " << cellpos(0) << "\n"
410
+				<< "Region limits: " << this->Origin(0) << " to " <<
411
+				this->Nx3(0)*this->Dx3(0)+this->Origin(0) << std::endl;
412
+			exit(EXIT_FAILURE);
413
+		}
414
+		
415
+		if (cellpos(1) < this->Origin(1) || cellpos(1) >= 
416
+			this->Nx3(1)*this->Dx3(1) + this->Origin(1)) {
417
+			std::cout << "Queried point exists outside model region.\n";
418
+			std::cout << "Given location: " << cellpos(1) << "\n"
419
+				<< "Region limits: " << this->Origin(1) << " to " <<
420
+				this->Nx3(1)*this->Dx3(1)+this->Origin(1) << std::endl;
421
+			exit(EXIT_FAILURE);
422
+		}
423
+
424
+		if (cellpos(2) < this->Origin(2) || cellpos(2) >= 
425
+			this->Nx3(2)*this->Dx3(2) + this->Origin(2)) {
426
+			std::cout << "Queried point exists outside model region.\n";
427
+			std::cout << "Given location: " << cellpos(2) << "\n"
428
+				<< "Region limits: " << this->Origin(2) << " to " <<
429
+				this->Nx3(2)*this->Dx3(2)+this->Origin(2) << std::endl;
430
+			exit(EXIT_FAILURE);
431
+		}
432
+
433
+		//Subtract off the origin just to make the maths easier
434
+		for (int ii=0;ii<3;ii++) {
435
+			cellposmo(ii) = cellpos(ii) - this->Origin(ii);
436
+		}
437
+
438
+		vecpos = 42;
439
+		/// \todo Finish this up.  Returning 42 right now.
440
+		return vecpos;		
441
+	}
442
+
443
+	Eigen::Matrix<int,3,1> StructuredGrid::GetNCells() {
444
+		return this->Nx3;
445
+	}
446
+
447
+	Vector3r StructuredGrid::GetOrigin() {
448
+		return this->Origin;
449
+	}
450
+
451
+	Vector3Xr StructuredGrid::GetCentroids() {
452
+		return this->VecCentroid;
453
+	}
454
+
455
+	Vector3r StructuredGrid::GetCentroid(const int& vecpos) {
456
+		return this->VecCentroid.row(vecpos);
457
+	}
458
+
459
+
460
+
461
+
462
+
463
+
464
+
465
+
466
+
467
+
468
+
469
+
470
+}

Loading…
Cancel
Save