/[escript]/trunk/dataexporter/src/FinleyMesh.cpp
ViewVC logotype

Contents of /trunk/dataexporter/src/FinleyMesh.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2810 - (show annotations)
Mon Dec 7 04:13:49 2009 UTC (10 years, 11 months ago) by caltinay
File size: 9011 byte(s)
Reincarnation of the escriptreader as a more flexible escriptexport library:
 - can be initialised with instances of escript::Data and Finley_Mesh
 - is now MPI aware at the EscriptDataset level including Silo writer
 - now uses boost shared pointers
The lib is currently only used by the escriptconvert tool but this is going to
change...

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14 #include <escriptexport/FinleyMesh.h>
15 #include <escriptexport/ElementData.h>
16 #include <escriptexport/NodeData.h>
17
18 #include <finley/CppAdapter/MeshAdapter.h>
19 extern "C" {
20 #include <finley/Mesh.h>
21 }
22
23 #include <iostream>
24
25 #if USE_NETCDF
26 #include <netcdf.hh>
27 #endif
28
29 #if USE_SILO
30 #include <silo.h>
31 #endif
32
33 using namespace std;
34
35 namespace escriptexport {
36
37 //
38 //
39 //
40 FinleyMesh::FinleyMesh() :
41 initialized(false)
42 {
43 }
44
45 //
46 //
47 //
48 FinleyMesh::FinleyMesh(const FinleyMesh& m)
49 {
50 nodes = NodeData_ptr(new NodeData(*m.nodes));
51 cells = ElementData_ptr(new ElementData(*m.cells));
52 faces = ElementData_ptr(new ElementData(*m.faces));
53 contacts = ElementData_ptr(new ElementData(*m.contacts));
54 initialized = m.initialized;
55 }
56
57 //
58 //
59 //
60 FinleyMesh::~FinleyMesh()
61 {
62 cleanup();
63 }
64
65 //
66 //
67 //
68 void FinleyMesh::cleanup()
69 {
70 nodes.reset();
71 cells.reset();
72 faces.reset();
73 contacts.reset();
74 initialized = false;
75 }
76
77 //
78 //
79 //
80 bool FinleyMesh::initFromEscript(escript::const_Domain_ptr escriptDomain)
81 {
82 cleanup();
83
84 finleyMesh = dynamic_cast<const finley::MeshAdapter*>(escriptDomain.get())->getFinley_Mesh();
85 if (!finleyMesh) {
86 return false;
87 }
88
89 nodes = NodeData_ptr(new NodeData("Elements"));
90 cells = ElementData_ptr(new ElementData("Elements", nodes));
91 faces = ElementData_ptr(new ElementData("FaceElements", nodes));
92 contacts = ElementData_ptr(new ElementData("ContactElements", nodes));
93
94 if (nodes->initFromFinley(finleyMesh->Nodes) &&
95 cells->initFromFinley(finleyMesh->Elements) &&
96 faces->initFromFinley(finleyMesh->FaceElements) &&
97 contacts->initFromFinley(finleyMesh->ContactElements)) {
98 initialized = true;
99 }
100
101 return initialized;
102 }
103
104 //
105 // Reads mesh and element data from NetCDF file with given name
106 //
107 bool FinleyMesh::initFromNetCDF(const string& filename)
108 {
109 cleanup();
110
111 #if USE_NETCDF
112 NcError ncerr(NcError::silent_nonfatal);
113 NcFile* input;
114
115 input = new NcFile(filename.c_str());
116 if (!input->is_valid()) {
117 cerr << "Could not open input file " << filename << "." << endl;
118 delete input;
119 return false;
120 }
121
122 nodes = NodeData_ptr(new NodeData("Elements"));
123 if (!nodes->readFromNc(input))
124 return false;
125
126 // Read all element types
127 cells = ElementData_ptr(new ElementData("Elements", nodes));
128 cells->readFromNc(input);
129 faces = ElementData_ptr(new ElementData("FaceElements", nodes));
130 faces->readFromNc(input);
131 contacts = ElementData_ptr(new ElementData("ContactElements", nodes));
132 contacts->readFromNc(input);
133
134 delete input;
135 initialized = true;
136 #endif
137
138 return initialized;
139 }
140
141 //
142 //
143 //
144 NodeData_ptr FinleyMesh::getMeshForFinleyFS(int functionSpace) const
145 {
146 NodeData_ptr result;
147
148 if (!initialized)
149 return result;
150
151 ElementData_ptr elements = getElementsForFinleyFS(functionSpace);
152 if (elements != NULL) {
153 if (functionSpace == FINLEY_NODES)
154 result = elements->getNodeMesh();
155 else {
156 if (elements->getReducedNumElements() > 0)
157 result = elements->getReducedNodeMesh();
158 else
159 result = elements->getNodeMesh();
160 }
161 }
162
163 return result;
164 }
165
166 //
167 //
168 //
169 ElementData_ptr FinleyMesh::getElementsForFinleyFS(int functionSpace) const
170 {
171 ElementData_ptr result;
172
173 if (!initialized) {
174 return result;
175 }
176
177 switch (functionSpace) {
178 case FINLEY_NODES:
179 case FINLEY_REDUCED_NODES:
180 case FINLEY_REDUCED_ELEMENTS:
181 case FINLEY_ELEMENTS:
182 result = cells;
183 break;
184
185 case FINLEY_REDUCED_FACE_ELEMENTS:
186 case FINLEY_FACE_ELEMENTS:
187 result = faces;
188 break;
189
190 case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
191 case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
192 case FINLEY_CONTACT_ELEMENTS_1:
193 case FINLEY_CONTACT_ELEMENTS_2:
194 result = contacts;
195 break;
196
197 default:
198 cerr << "Unsupported function space type " << functionSpace
199 << "!" << endl;
200 }
201 return result;
202 }
203
204 //
205 // Returns a vector of strings containing Silo mesh names that have been
206 // written
207 //
208 StringVec FinleyMesh::getMeshNames() const
209 {
210 StringVec res;
211 if (initialized) {
212 StringVec tmpVec;
213 tmpVec = cells->getMeshNames();
214 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
215 tmpVec = faces->getMeshNames();
216 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
217 tmpVec = contacts->getMeshNames();
218 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
219 }
220 return res;
221 }
222
223 //
224 //
225 //
226 NodeData_ptr FinleyMesh::getMeshByName(const string& name) const
227 {
228 NodeData_ptr ret;
229 if (!initialized) {
230 return ret;
231 }
232 if (name == "Elements")
233 ret = cells->getNodeMesh();
234 else if (name == "ReducedElements")
235 ret = cells->getReducedNodeMesh();
236 else if (name == "FaceElements")
237 ret = faces->getNodeMesh();
238 else if (name == "ReducedFaceElements")
239 ret = faces->getReducedNodeMesh();
240 else if (name == "ContactElements")
241 ret = contacts->getNodeMesh();
242 else if (name == "ReducedContactElements")
243 ret = contacts->getReducedNodeMesh();
244
245 return ret;
246 }
247
248 //
249 // Returns a vector of strings containing Silo variable names that have
250 // been written
251 //
252 StringVec FinleyMesh::getVarNames() const
253 {
254 StringVec res;
255
256 if (initialized) {
257 res = nodes->getVarNames();
258 StringVec tmpVec = cells->getVarNames();
259 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
260 tmpVec = faces->getVarNames();
261 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
262 tmpVec = contacts->getVarNames();
263 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
264 }
265
266 return res;
267 }
268
269 //
270 //
271 //
272 const IntVec& FinleyMesh::getVarDataByName(const string& name) const
273 {
274 if (!initialized) {
275 throw "Mesh not initialized";
276 }
277
278 if (name.compare(0, 6, "Nodes_") == 0)
279 return nodes->getVarDataByName(name);
280 else if (name.compare(0, 9, "Elements_") == 0)
281 return cells->getVarDataByName(name);
282 else if (name.compare(0, 13, "FaceElements_") == 0)
283 return faces->getVarDataByName(name);
284 else if (name.compare(0, 16, "ContactElements_") == 0)
285 return contacts->getVarDataByName(name);
286 else
287 throw "Invalid variable name";
288 }
289
290 //
291 //
292 //
293 ElementData_ptr FinleyMesh::getElementsByName(const string& name) const
294 {
295 ElementData_ptr ret;
296 if (name == "Elements" || name == "ReducedElements")
297 ret = cells;
298 else if (name == "FaceElements" || name == "ReducedFaceElements")
299 ret = faces;
300 else if (name == "ContactElements" || name == "ReducedContactElements")
301 ret = contacts;
302
303 return ret;
304 }
305
306 //
307 //
308 //
309 void FinleyMesh::reorderGhostZones(int ownIndex)
310 {
311 if (initialized) {
312 cells->reorderGhostZones(ownIndex);
313 faces->reorderGhostZones(ownIndex);
314 contacts->reorderGhostZones(ownIndex);
315 #ifdef _DEBUG
316 cout << "block " << ownIndex << " has " << cells->getGhostCount()
317 << " (reduced=" << cells->getReducedGhostCount()
318 << ") ghost zones," << endl;
319 cout << "\t" << faces->getGhostCount() << " ("
320 << faces->getReducedGhostCount() << ") ghost faces," << endl;
321 cout << "\t" << contacts->getGhostCount() << " ("
322 << contacts->getReducedGhostCount() << ") ghost contacts." << endl;
323 #endif
324 }
325 }
326
327 //
328 //
329 //
330 void FinleyMesh::removeGhostZones(int ownIndex)
331 {
332 if (initialized) {
333 cells->removeGhostZones(ownIndex);
334 faces->removeGhostZones(ownIndex);
335 contacts->removeGhostZones(ownIndex);
336 #ifdef _DEBUG
337 cout << "After removing ghost zones there are" << endl;
338 cout << " " << nodes->getNumNodes() << " Nodes, ";
339 cout << cells->getCount() << " Elements, ";
340 cout << faces->getCount() << " Face elements, ";
341 cout << contacts->getCount() << " Contact elements left." << endl;
342 #endif
343 }
344 }
345
346 //
347 //
348 //
349 bool FinleyMesh::writeToSilo(DBfile* dbfile, const string& pathInSilo)
350 {
351 #if USE_SILO
352 // Write nodes, elements and mesh variables
353 if (!initialized ||
354 !cells->writeToSilo(dbfile, pathInSilo) ||
355 !faces->writeToSilo(dbfile, pathInSilo) ||
356 !contacts->writeToSilo(dbfile, pathInSilo))
357 return false;
358
359 siloPath = pathInSilo;
360 return true;
361
362 #else // !USE_SILO
363 return false;
364 #endif
365 }
366
367 } // namespace escriptexport
368

  ViewVC Help
Powered by ViewVC 1.1.26