/[escript]/branches/split/weipa/src/FinleyDomain.cpp
ViewVC logotype

Contents of /branches/split/weipa/src/FinleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4724 - (show annotations)
Thu Mar 6 05:22:12 2014 UTC (5 years, 1 month ago) by jfenwick
File size: 12085 byte(s)
Work towards parallel domains

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

  ViewVC Help
Powered by ViewVC 1.1.26