/[escript]/trunk/weipa/src/FinleyDomain.cpp
ViewVC logotype

Contents of /trunk/weipa/src/FinleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4492 - (show annotations)
Tue Jul 2 01:44:11 2013 UTC (5 years, 9 months ago) by caltinay
File size: 12016 byte(s)
Finley changes that were held back while in release mode - moved more stuff
into finley namespace.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16 #include <weipa/FinleyDomain.h>
17 #include <weipa/FinleyNodes.h>
18 #include <weipa/DataVar.h>
19
20 #ifndef VISIT_PLUGIN
21 #include <dudley/CppAdapter/MeshAdapter.h>
22 #include <finley/CppAdapter/MeshAdapter.h>
23 #include <dudley/Mesh.h>
24 #include <finley/Mesh.h>
25 #endif
26
27 #include <iostream>
28
29 #if USE_NETCDF
30 #include <netcdfcpp.h>
31 #endif
32
33 #if USE_SILO
34 #include <silo.h>
35 #endif
36
37 using namespace std;
38
39 namespace weipa {
40
41 //
42 //
43 //
44 FinleyDomain::FinleyDomain() :
45 initialized(false)
46 {
47 }
48
49 //
50 //
51 //
52 FinleyDomain::FinleyDomain(const FinleyDomain& m) :
53 boost::enable_shared_from_this<FinleyDomain>()
54 {
55 nodes = FinleyNodes_ptr(new FinleyNodes(*m.nodes));
56 cells = FinleyElements_ptr(new FinleyElements(*m.cells));
57 faces = FinleyElements_ptr(new FinleyElements(*m.faces));
58 contacts = FinleyElements_ptr(new FinleyElements(*m.contacts));
59 initialized = m.initialized;
60 }
61
62 //
63 //
64 //
65 FinleyDomain::~FinleyDomain()
66 {
67 cleanup();
68 }
69
70 //
71 //
72 //
73 void FinleyDomain::cleanup()
74 {
75 nodes.reset();
76 cells.reset();
77 faces.reset();
78 contacts.reset();
79 initialized = false;
80 }
81
82 //
83 //
84 //
85 bool FinleyDomain::initFromEscript(const escript::AbstractDomain* escriptDomain)
86 {
87 #ifndef VISIT_PLUGIN
88 cleanup();
89
90 if (dynamic_cast<const finley::MeshAdapter*>(escriptDomain)) {
91 const Finley_Mesh* finleyMesh =
92 dynamic_cast<const finley::MeshAdapter*>(escriptDomain)
93 ->getFinley_Mesh();
94
95 nodes = FinleyNodes_ptr(new FinleyNodes("Elements"));
96 cells = FinleyElements_ptr(new FinleyElements("Elements", nodes));
97 faces = FinleyElements_ptr(new FinleyElements("FaceElements", nodes));
98 contacts = FinleyElements_ptr(new FinleyElements("ContactElements", nodes));
99
100 if (nodes->initFromFinley(finleyMesh->Nodes) &&
101 cells->initFromFinley(finleyMesh->Elements) &&
102 faces->initFromFinley(finleyMesh->FaceElements) &&
103 contacts->initFromFinley(finleyMesh->ContactElements)) {
104 initialized = true;
105 }
106 } else if (dynamic_cast<const dudley::MeshAdapter*>(escriptDomain)) {
107 const Dudley_Mesh* dudleyMesh =
108 dynamic_cast<const dudley::MeshAdapter*>(escriptDomain)
109 ->getDudley_Mesh();
110
111 nodes = FinleyNodes_ptr(new FinleyNodes("Elements"));
112 cells = FinleyElements_ptr(new FinleyElements("Elements", nodes));
113 faces = FinleyElements_ptr(new FinleyElements("FaceElements", nodes));
114 contacts = FinleyElements_ptr(new FinleyElements("ContactElements", nodes));
115
116 if (nodes->initFromDudley(dudleyMesh->Nodes) &&
117 cells->initFromDudley(dudleyMesh->Elements) &&
118 faces->initFromDudley(dudleyMesh->FaceElements)) {
119 initialized = true;
120 }
121 }
122
123 return initialized;
124 #else // VISIT_PLUGIN
125 return false;
126 #endif
127 }
128
129 //
130 // Reads mesh and element data from NetCDF file with given name
131 //
132 bool FinleyDomain::initFromFile(const string& filename)
133 {
134 cleanup();
135
136 #if USE_NETCDF
137 NcError ncerr(NcError::silent_nonfatal);
138 NcFile* input;
139
140 input = new NcFile(filename.c_str());
141 if (!input->is_valid()) {
142 cerr << "Could not open input file " << filename << "." << endl;
143 delete input;
144 return false;
145 }
146
147 nodes = FinleyNodes_ptr(new FinleyNodes("Elements"));
148 if (!nodes->readFromNc(input))
149 return false;
150
151 // Read all element types
152 cells = FinleyElements_ptr(new FinleyElements("Elements", nodes));
153 cells->readFromNc(input);
154 faces = FinleyElements_ptr(new FinleyElements("FaceElements", nodes));
155 faces->readFromNc(input);
156 contacts = FinleyElements_ptr(new FinleyElements("ContactElements", nodes));
157 contacts->readFromNc(input);
158
159 delete input;
160 initialized = true;
161 #endif
162
163 return initialized;
164 }
165
166 Centering FinleyDomain::getCenteringForFunctionSpace(int fsCode) const
167 {
168 return (fsCode==FINLEY_REDUCED_NODES || fsCode==FINLEY_NODES ?
169 NODE_CENTERED : ZONE_CENTERED);
170 }
171
172 //
173 //
174 //
175 NodeData_ptr FinleyDomain::getMeshForFunctionSpace(int fsCode) const
176 {
177 NodeData_ptr result;
178
179 if (!initialized)
180 return result;
181
182 ElementData_ptr elements = getElementsForFunctionSpace(fsCode);
183 if (elements != NULL)
184 result = elements->getNodes();
185
186 return result;
187 }
188
189 //
190 //
191 //
192 ElementData_ptr FinleyDomain::getElementsForFunctionSpace(int fsCode) const
193 {
194 ElementData_ptr result;
195
196 if (!initialized) {
197 return result;
198 }
199
200 if (fsCode == FINLEY_NODES) {
201 result = cells;
202 } else if (fsCode == FINLEY_REDUCED_NODES) {
203 result = cells->getReducedElements();
204 if (!result)
205 result = cells;
206 } else {
207 switch (fsCode) {
208 case FINLEY_REDUCED_ELEMENTS:
209 case FINLEY_ELEMENTS:
210 result = cells;
211 break;
212
213 case FINLEY_REDUCED_FACE_ELEMENTS:
214 case FINLEY_FACE_ELEMENTS:
215 result = faces;
216 break;
217
218 case FINLEY_REDUCED_CONTACT_ELEMENTS_1:
219 case FINLEY_REDUCED_CONTACT_ELEMENTS_2:
220 case FINLEY_CONTACT_ELEMENTS_1:
221 case FINLEY_CONTACT_ELEMENTS_2:
222 result = contacts;
223 break;
224
225 default: {
226 cerr << "Unsupported function space type " << fsCode
227 << "!" << endl;
228 return result;
229 }
230 }
231 int typeId = static_cast<FinleyElements*>(result.get())
232 ->getFinleyTypeId();
233 if (typeId != finley::Line3Macro &&
234 typeId != finley::Rec9 && typeId != finley::Rec9Macro &&
235 typeId != finley::Hex27 &&
236 typeId != finley::Hex27Macro && typeId != finley::Tri6 &&
237 typeId != finley::Tri6Macro && typeId != finley::Tet10 &&
238 typeId != finley::Tet10Macro) {
239 if (result->getReducedElements())
240 result = result->getReducedElements();
241 }
242 }
243
244 return result;
245 }
246
247 //
248 // Returns a vector of strings containing mesh names for this domain
249 //
250 StringVec FinleyDomain::getMeshNames() const
251 {
252 StringVec res;
253 if (initialized) {
254 StringVec tmpVec;
255 tmpVec = cells->getMeshNames();
256 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
257 tmpVec = faces->getMeshNames();
258 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
259 tmpVec = contacts->getMeshNames();
260 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
261 }
262 return res;
263 }
264
265 //
266 // Returns a vector of strings containing mesh variable names for this domain
267 //
268 StringVec FinleyDomain::getVarNames() const
269 {
270 StringVec res;
271
272 if (initialized) {
273 res = nodes->getVarNames();
274 StringVec tmpVec = cells->getVarNames();
275 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
276 tmpVec = faces->getVarNames();
277 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
278 tmpVec = contacts->getVarNames();
279 res.insert(res.end(), tmpVec.begin(), tmpVec.end());
280 }
281
282 return res;
283 }
284
285 //
286 //
287 //
288 DataVar_ptr FinleyDomain::getDataVarByName(const string& name) const
289 {
290 if (!initialized) {
291 throw "Domain not initialized";
292 }
293
294 DataVar_ptr var(new DataVar(name));
295 if (name.find("ContactElements_") != name.npos) {
296 const IntVec& data = contacts->getVarDataByName(name);
297 string elementName = name.substr(0, name.find('_'));
298 ElementData_ptr elements = getElementsByName(elementName);
299 var->initFromMeshData(shared_from_this(), data,
300 FINLEY_CONTACT_ELEMENTS_1, ZONE_CENTERED, elements->getNodes(),
301 elements->getIDs());
302 } else if (name.find("FaceElements_") != name.npos) {
303 const IntVec& data = faces->getVarDataByName(name);
304 string elementName = name.substr(0, name.find('_'));
305 ElementData_ptr elements = getElementsByName(elementName);
306 var->initFromMeshData(shared_from_this(), data,
307 FINLEY_FACE_ELEMENTS, ZONE_CENTERED, elements->getNodes(),
308 elements->getIDs());
309 } else if (name.find("Elements_") != name.npos) {
310 const IntVec& data = cells->getVarDataByName(name);
311 string elementName = name.substr(0, name.find('_'));
312 ElementData_ptr elements = getElementsByName(elementName);
313 var->initFromMeshData(shared_from_this(), data, FINLEY_ELEMENTS,
314 ZONE_CENTERED, elements->getNodes(), elements->getIDs());
315 } else if (name.find("Nodes_") != name.npos) {
316 const IntVec& data = nodes->getVarDataByName(name);
317 var->initFromMeshData(shared_from_this(), data, FINLEY_NODES,
318 NODE_CENTERED, getNodes(), getNodes()->getNodeIDs());
319 } else {
320 cerr << "WARNING: Unrecognized domain variable '" << name << "'\n";
321 return DataVar_ptr();
322 }
323
324 return var;
325 }
326
327 //
328 //
329 //
330 ElementData_ptr FinleyDomain::getElementsByName(const string& name) const
331 {
332 ElementData_ptr ret;
333 if (name == "Elements")
334 ret = cells;
335 else if (name == "ReducedElements")
336 ret = cells->getReducedElements();
337 else if (name == "FaceElements")
338 ret = faces;
339 else if (name == "ReducedFaceElements")
340 ret = faces->getReducedElements();
341 else if (name == "ContactElements")
342 ret = contacts;
343 else if (name == "ReducedContactElements")
344 ret = contacts->getReducedElements();
345
346 return ret;
347 }
348
349 //
350 //
351 //
352 NodeData_ptr FinleyDomain::getMeshByName(const string& name) const
353 {
354 NodeData_ptr ret;
355 if (initialized) {
356 ElementData_ptr els = getElementsByName(name);
357 if (els)
358 ret = els->getNodes();
359 }
360
361 return ret;
362 }
363
364 //
365 //
366 //
367 void FinleyDomain::reorderGhostZones(int ownIndex)
368 {
369 if (initialized) {
370 cells->reorderGhostZones(ownIndex);
371 faces->reorderGhostZones(ownIndex);
372 contacts->reorderGhostZones(ownIndex);
373 #ifdef _DEBUG
374 cout << "block " << ownIndex << " has " << cells->getGhostCount()
375 << " ghost zones," << endl;
376 cout << "\t" << faces->getGhostCount() << " ghost faces," << endl;
377 cout << "\t" << contacts->getGhostCount() << " ghost contacts." << endl;
378 #endif
379 }
380 }
381
382 //
383 //
384 //
385 void FinleyDomain::removeGhostZones(int ownIndex)
386 {
387 if (initialized) {
388 cells->removeGhostZones(ownIndex);
389 faces->removeGhostZones(ownIndex);
390 contacts->removeGhostZones(ownIndex);
391 #ifdef _DEBUG
392 cout << "After removing ghost zones there are" << endl;
393 cout << " " << nodes->getNumNodes() << " Nodes, ";
394 cout << cells->getCount() << " Elements, ";
395 cout << faces->getCount() << " Face elements, ";
396 cout << contacts->getCount() << " Contact elements left." << endl;
397 #endif
398 }
399 }
400
401 //
402 //
403 //
404 bool FinleyDomain::writeToSilo(DBfile* dbfile, const string& pathInSilo,
405 const StringVec& labels, const StringVec& units,
406 bool writeMeshData)
407 {
408 #if USE_SILO
409 // Write nodes, elements and mesh variables
410 if (!initialized ||
411 !cells->writeToSilo(dbfile, pathInSilo, labels, units, writeMeshData) ||
412 !faces->writeToSilo(dbfile, pathInSilo, labels, units, writeMeshData) ||
413 !contacts->writeToSilo(dbfile, pathInSilo, labels, units, writeMeshData))
414 return false;
415
416 siloPath = pathInSilo;
417 return true;
418
419 #else // !USE_SILO
420 return false;
421 #endif
422 }
423
424 } // namespace weipa
425

  ViewVC Help
Powered by ViewVC 1.1.26