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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2810 - (hide 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 caltinay 2183
2     /*******************************************************
3     *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 caltinay 2183 * 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 caltinay 2810 #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 caltinay 2183 #include <netcdf.hh>
27 caltinay 2810 #endif
28    
29     #if USE_SILO
30 caltinay 2183 #include <silo.h>
31     #endif
32    
33     using namespace std;
34    
35 caltinay 2810 namespace escriptexport {
36 caltinay 2187
37 caltinay 2183 //
38     //
39     //
40 caltinay 2810 FinleyMesh::FinleyMesh() :
41     initialized(false)
42 caltinay 2183 {
43     }
44    
45     //
46     //
47     //
48 caltinay 2810 FinleyMesh::FinleyMesh(const FinleyMesh& m)
49 caltinay 2183 {
50 caltinay 2810 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 caltinay 2183 }
56    
57     //
58     //
59     //
60 caltinay 2810 FinleyMesh::~FinleyMesh()
61 caltinay 2183 {
62 caltinay 2810 cleanup();
63 caltinay 2183 }
64    
65     //
66 caltinay 2810 //
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 caltinay 2183 // Returns a vector of strings containing Silo mesh names that have been
206     // written
207     //
208 caltinay 2810 StringVec FinleyMesh::getMeshNames() const
209 caltinay 2183 {
210     StringVec res;
211 caltinay 2810 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 caltinay 2183 return res;
221     }
222    
223 caltinay 2810 //
224     //
225     //
226     NodeData_ptr FinleyMesh::getMeshByName(const string& name) const
227 caltinay 2183 {
228 caltinay 2810 NodeData_ptr ret;
229     if (!initialized) {
230     return ret;
231     }
232 caltinay 2183 if (name == "Elements")
233 caltinay 2810 ret = cells->getNodeMesh();
234 caltinay 2183 else if (name == "ReducedElements")
235 caltinay 2810 ret = cells->getReducedNodeMesh();
236 caltinay 2183 else if (name == "FaceElements")
237 caltinay 2810 ret = faces->getNodeMesh();
238 caltinay 2183 else if (name == "ReducedFaceElements")
239 caltinay 2810 ret = faces->getReducedNodeMesh();
240 caltinay 2183 else if (name == "ContactElements")
241 caltinay 2810 ret = contacts->getNodeMesh();
242 caltinay 2183 else if (name == "ReducedContactElements")
243 caltinay 2810 ret = contacts->getReducedNodeMesh();
244 caltinay 2183
245     return ret;
246     }
247    
248     //
249     // Returns a vector of strings containing Silo variable names that have
250     // been written
251     //
252 caltinay 2810 StringVec FinleyMesh::getVarNames() const
253 caltinay 2183 {
254     StringVec res;
255 caltinay 2810
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 caltinay 2183
266     return res;
267     }
268    
269 caltinay 2810 //
270     //
271     //
272     const IntVec& FinleyMesh::getVarDataByName(const string& name) const
273 caltinay 2183 {
274 caltinay 2810 if (!initialized) {
275     throw "Mesh not initialized";
276     }
277    
278     if (name.compare(0, 6, "Nodes_") == 0)
279     return nodes->getVarDataByName(name);
280 caltinay 2183 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 caltinay 2201 throw "Invalid variable name";
288 caltinay 2183 }
289    
290 caltinay 2810 //
291     //
292     //
293     ElementData_ptr FinleyMesh::getElementsByName(const string& name) const
294 caltinay 2183 {
295 caltinay 2810 ElementData_ptr ret;
296 caltinay 2183 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 caltinay 2810 void FinleyMesh::reorderGhostZones(int ownIndex)
310 caltinay 2183 {
311 caltinay 2810 if (initialized) {
312     cells->reorderGhostZones(ownIndex);
313     faces->reorderGhostZones(ownIndex);
314     contacts->reorderGhostZones(ownIndex);
315 caltinay 2183 #ifdef _DEBUG
316 caltinay 2810 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 caltinay 2183 #endif
324 caltinay 2810 }
325 caltinay 2183 }
326    
327     //
328     //
329     //
330 caltinay 2810 void FinleyMesh::removeGhostZones(int ownIndex)
331 caltinay 2183 {
332 caltinay 2810 if (initialized) {
333     cells->removeGhostZones(ownIndex);
334     faces->removeGhostZones(ownIndex);
335     contacts->removeGhostZones(ownIndex);
336 caltinay 2183 #ifdef _DEBUG
337 caltinay 2810 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 caltinay 2183 #endif
343 caltinay 2810 }
344 caltinay 2183 }
345    
346     //
347     //
348     //
349 caltinay 2810 bool FinleyMesh::writeToSilo(DBfile* dbfile, const string& pathInSilo)
350 caltinay 2183 {
351 caltinay 2810 #if USE_SILO
352     // Write nodes, elements and mesh variables
353     if (!initialized ||
354     !cells->writeToSilo(dbfile, pathInSilo) ||
355 caltinay 2183 !faces->writeToSilo(dbfile, pathInSilo) ||
356 caltinay 2810 !contacts->writeToSilo(dbfile, pathInSilo))
357 caltinay 2183 return false;
358    
359 caltinay 2810 siloPath = pathInSilo;
360     return true;
361 caltinay 2183
362 caltinay 2810 #else // !USE_SILO
363 caltinay 2183 return false;
364     #endif
365     }
366    
367 caltinay 2810 } // namespace escriptexport
368 caltinay 2187

  ViewVC Help
Powered by ViewVC 1.1.26