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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3247 - (show annotations)
Wed Oct 6 05:53:06 2010 UTC (8 years, 6 months ago) by caltinay
File size: 10831 byte(s)
Fixed name clashes between dudley and finley so both can be used
simultaneously.

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

  ViewVC Help
Powered by ViewVC 1.1.26