/[escript]/trunk/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Contents of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1345 - (show annotations)
Wed Nov 14 07:53:34 2007 UTC (11 years, 7 months ago) by ksteube
File size: 80145 byte(s)
Created LoadMesh to read a mesh from a distributed NetCDF file.
Can read nodes but not elements yet.

1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 #include "MeshAdapter.h"
17 #include "escript/Data.h"
18 #include "escript/DataFactory.h"
19 #ifdef USE_NETCDF
20 #include <netcdfcpp.h>
21 #endif
22 extern "C" {
23 #include "escript/blocktimer.h"
24 }
25 #include <vector>
26
27 using namespace std;
28 using namespace escript;
29
30 namespace finley {
31
32 //
33 // define the static constants
34 MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
35 const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
36 const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
37 const int MeshAdapter::Nodes=FINLEY_NODES;
38 const int MeshAdapter::ReducedNodes=FINLEY_REDUCED_NODES;
39 const int MeshAdapter::Elements=FINLEY_ELEMENTS;
40 const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;
41 const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
42 const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;
43 const int MeshAdapter::Points=FINLEY_POINTS;
44 const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
45 const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;
46 const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
47 const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
48
49 MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
50 {
51 setFunctionSpaceTypeNames();
52 //
53 // need to use a null_deleter as Finley_Mesh_free deletes the pointer
54 // for us.
55 m_finleyMesh.reset(finleyMesh,null_deleter());
56 }
57
58 //
59 // The copy constructor should just increment the use count
60 MeshAdapter::MeshAdapter(const MeshAdapter& in):
61 m_finleyMesh(in.m_finleyMesh)
62 {
63 setFunctionSpaceTypeNames();
64 }
65
66 MeshAdapter::~MeshAdapter()
67 {
68 //
69 // I hope the case for the pointer being zero has been taken care of.
70 // cout << "In MeshAdapter destructor." << endl;
71 if (m_finleyMesh.unique()) {
72 Finley_Mesh_free(m_finleyMesh.get());
73 }
74 }
75
76 int MeshAdapter::getMPISize() const
77 {
78 return m_finleyMesh.get()->MPIInfo->size;
79 }
80 int MeshAdapter::getMPIRank() const
81 {
82 return m_finleyMesh.get()->MPIInfo->rank;
83 }
84
85
86 Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
87 return m_finleyMesh.get();
88 }
89
90 void MeshAdapter::write(const std::string& fileName) const
91 {
92 char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
93 strcpy(fName,fileName.c_str());
94 Finley_Mesh_write(m_finleyMesh.get(),fName);
95 checkFinleyError();
96 TMPMEMFREE(fName);
97 }
98
99 void MeshAdapter::Print_Mesh_Info(const bool full=false) const
100 {
101 Finley_PrintMesh_Info(m_finleyMesh.get(), full);
102 }
103
104 void MeshAdapter::dump(const std::string& fileName) const
105 {
106 #ifdef USE_NETCDF
107 const NcDim* ncdims[25]; /* ksteube how big should ncdims be? */
108 NcVar *ids, *data;
109 int *int_ptr;
110 Finley_Mesh *mesh = m_finleyMesh.get();
111 int mpi_size = mesh->MPIInfo->size;
112 int mpi_rank = mesh->MPIInfo->rank;
113 int numDim = mesh->Nodes->numDim;
114 int numNodes = mesh->Nodes->numNodes;
115 int num_Elements = mesh->Elements->numElements;
116 int num_FaceElements = mesh->FaceElements->numElements;
117 int num_ContactElements = mesh->ContactElements->numElements;
118 int num_Points = mesh->Points->numElements;
119 int num_Elements_numNodes = mesh->Elements->numNodes;
120 int num_FaceElements_numNodes = mesh->FaceElements->numNodes;
121 int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
122 char *newFileName = Paso_MPI_appendRankToFileName(strdup(fileName.c_str()), mpi_size, mpi_rank);
123
124 // NetCDF error handler
125 NcError err(NcError::verbose_nonfatal);
126 // Create the file.
127 NcFile dataFile(newFileName, NcFile::Replace);
128 // check if writing was successful
129 if (!dataFile.is_valid())
130 throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);
131
132 // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)
133 if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
134 throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);
135 if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
136 throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);
137 if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
138 throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);
139 if (num_Elements>0)
140 if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
141 throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);
142 if (num_FaceElements>0)
143 if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
144 throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);
145 if (num_ContactElements>0)
146 if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
147 throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);
148 if (num_Points>0)
149 if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
150 throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);
151 if (num_Elements>0)
152 if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
153 throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);
154 if (num_FaceElements>0)
155 if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
156 throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);
157 if (num_ContactElements>0)
158 if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
159 throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);
160
161 // Attributes: MPI size, MPI rank, Name, order, reduced_order
162 if (!dataFile.add_att("mpi_size", mpi_size) )
163 throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);
164 if (!dataFile.add_att("mpi_rank", mpi_rank) )
165 throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);
166 if (!dataFile.add_att("Name",mesh->Name) )
167 throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);
168 if (!dataFile.add_att("numDim",numDim) )
169 throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);
170 if (!dataFile.add_att("order",mesh->order) )
171 throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);
172 if (!dataFile.add_att("reduced_order",mesh->reduced_order) )
173 throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);
174 if (!dataFile.add_att("numNodes",numNodes) )
175 throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);
176 if (!dataFile.add_att("num_Elements",num_Elements) )
177 throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);
178 if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
179 throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);
180 if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
181 throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);
182 if (!dataFile.add_att("num_Points",num_Points) )
183 throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);
184
185 // // // // // Nodes // // // // //
186
187 // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)
188 if (numNodes>0) {
189
190 // Nodes Id
191 if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
192 throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);
193 int_ptr = &mesh->Nodes->Id[0];
194 if (! (ids->put(int_ptr, numNodes)) )
195 throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);
196
197 // Nodes Tag
198 if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
199 throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);
200 int_ptr = &mesh->Nodes->Tag[0];
201 if (! (ids->put(int_ptr, numNodes)) )
202 throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);
203
204 // Nodes gDOF
205 if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
206 throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);
207 int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
208 if (! (ids->put(int_ptr, numNodes)) )
209 throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);
210
211 // Nodes global node index
212 if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
213 throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);
214 int_ptr = &mesh->Nodes->globalNodesIndex[0];
215 if (! (ids->put(int_ptr, numNodes)) )
216 throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);
217
218 // Nodes grDof
219 if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
220 throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);
221 int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
222 if (! (ids->put(int_ptr, numNodes)) )
223 throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);
224
225 // Nodes grNI
226 if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
227 throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);
228 int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
229 if (! (ids->put(int_ptr, numNodes)) )
230 throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);
231
232 // Nodes Coordinates
233 if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
234 throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);
235 if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
236 throw DataException("Error - MeshAdapter::dump: copy Nodes_Coordinates to netCDF buffer failed: " + *newFileName);
237
238 // Nodes degreesOfFreedomDistribution
239 if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
240 throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);
241 int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
242 if (! (ids->put(int_ptr, mpi_size+1)) )
243 throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);
244
245 }
246
247 // // // // // Elements // // // // //
248
249 if (num_Elements>0) {
250
251 // Temp storage to gather node IDs
252 int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
253
254 // Elements ReferenceElement->ElementTypeId
255 if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )
256 throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);
257
258 // Elements_Id
259 if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
260 throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);
261 int_ptr = &mesh->Elements->Id[0];
262 if (! (ids->put(int_ptr, num_Elements)) )
263 throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);
264
265 // Elements_Tag
266 if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
267 throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);
268 int_ptr = &mesh->Elements->Tag[0];
269 if (! (ids->put(int_ptr, num_Elements)) )
270 throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);
271
272 // Elements_Owner
273 if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
274 throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);
275 int_ptr = &mesh->Elements->Owner[0];
276 if (! (ids->put(int_ptr, num_Elements)) )
277 throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);
278
279 // Elements_Color
280 if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
281 throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);
282 int_ptr = &mesh->Elements->Color[0];
283 if (! (ids->put(int_ptr, num_Elements)) )
284 throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);
285
286 // Elements_Nodes
287 for (int i=0; i<num_Elements; i++)
288 for (int j=0; j<num_Elements_numNodes; j++)
289 Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)] = mesh->Nodes->Id[mesh->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]];
290 if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
291 throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);
292 if (! (ids->put(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes)) )
293 throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);
294
295 TMPMEMFREE(Elements_Nodes);
296
297 }
298
299 // // // // // Face_Elements // // // // //
300
301 if (num_FaceElements>0) {
302
303 // Temp storage to gather node IDs
304 int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
305
306 // FaceElements ReferenceElement->ElementTypeId
307 if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )
308 throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);
309
310 // FaceElements_Id
311 if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
312 throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);
313 int_ptr = &mesh->FaceElements->Id[0];
314 if (! (ids->put(int_ptr, num_FaceElements)) )
315 throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);
316
317 // FaceElements_Tag
318 if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
319 throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);
320 int_ptr = &mesh->FaceElements->Tag[0];
321 if (! (ids->put(int_ptr, num_FaceElements)) )
322 throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);
323
324 // FaceElements_Owner
325 if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
326 throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);
327 int_ptr = &mesh->FaceElements->Owner[0];
328 if (! (ids->put(int_ptr, num_FaceElements)) )
329 throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);
330
331 // FaceElements_Color
332 if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
333 throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);
334 int_ptr = &mesh->FaceElements->Color[0];
335 if (! (ids->put(int_ptr, num_FaceElements)) )
336 throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);
337
338 // FaceElements_Nodes
339 for (int i=0; i<num_FaceElements; i++)
340 for (int j=0; j<num_FaceElements_numNodes; j++)
341 FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = mesh->Nodes->Id[mesh->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]];
342 if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
343 throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);
344 if (! (ids->put(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
345 throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);
346
347 TMPMEMFREE(FaceElements_Nodes);
348
349 }
350
351 // // // // // Contact_Elements // // // // //
352
353 if (num_ContactElements>0) {
354
355 // Temp storage to gather node IDs
356 int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
357
358 // ContactElements ReferenceElement->ElementTypeId
359 if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )
360 throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);
361
362 // ContactElements_Id
363 if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
364 throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);
365 int_ptr = &mesh->ContactElements->Id[0];
366 if (! (ids->put(int_ptr, num_ContactElements)) )
367 throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);
368
369 // ContactElements_Tag
370 if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
371 throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);
372 int_ptr = &mesh->ContactElements->Tag[0];
373 if (! (ids->put(int_ptr, num_ContactElements)) )
374 throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);
375
376 // ContactElements_Owner
377 if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
378 throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);
379 int_ptr = &mesh->ContactElements->Owner[0];
380 if (! (ids->put(int_ptr, num_ContactElements)) )
381 throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);
382
383 // ContactElements_Color
384 if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
385 throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);
386 int_ptr = &mesh->ContactElements->Color[0];
387 if (! (ids->put(int_ptr, num_ContactElements)) )
388 throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);
389
390 // ContactElements_Nodes
391 for (int i=0; i<num_ContactElements; i++)
392 for (int j=0; j<num_ContactElements_numNodes; j++)
393 ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)] = mesh->Nodes->Id[mesh->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]];
394 if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
395 throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);
396 if (! (ids->put(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
397 throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);
398
399 TMPMEMFREE(ContactElements_Nodes);
400
401 }
402
403 // // // // // Points // // // // //
404
405 if (num_Points>0) {
406
407 fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");
408
409 // Temp storage to gather node IDs
410 int *Points_Nodes = TMPMEMALLOC(num_Points,int);
411
412 // Points ReferenceElement->ElementTypeId
413 if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )
414 throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);
415
416 // Points_Id
417 if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
418 throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);
419 int_ptr = &mesh->Points->Id[0];
420 if (! (ids->put(int_ptr, num_Points)) )
421 throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);
422
423 // Points_Tag
424 if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
425 throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);
426 int_ptr = &mesh->Points->Tag[0];
427 if (! (ids->put(int_ptr, num_Points)) )
428 throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);
429
430 // Points_Owner
431 if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
432 throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);
433 int_ptr = &mesh->Points->Owner[0];
434 if (! (ids->put(int_ptr, num_Points)) )
435 throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);
436
437 // Points_Color
438 if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
439 throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);
440 int_ptr = &mesh->Points->Color[0];
441 if (! (ids->put(int_ptr, num_Points)) )
442 throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);
443
444 // Points_Nodes
445 // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
446 for (int i=0; i<num_Points; i++)
447 Points_Nodes[i] = mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]];
448 if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
449 throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);
450 if (! (ids->put(&(Points_Nodes[0]), num_Points)) )
451 throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);
452
453 TMPMEMFREE(Points_Nodes);
454
455 }
456
457 // // // // // TagMap // // // // //
458
459 fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump does not save the TagMap yet\n\n\n");
460
461 // NetCDF file is closed by destructor of NcFile object
462 #else
463 Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
464 #endif /* USE_NETCDF */
465 checkFinleyError();
466 }
467
468 string MeshAdapter::getDescription() const
469 {
470 return "FinleyMesh";
471 }
472
473 string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
474 {
475 FunctionSpaceNamesMapType::iterator loc;
476 loc=m_functionSpaceTypeNames.find(functionSpaceType);
477 if (loc==m_functionSpaceTypeNames.end()) {
478 return "Invalid function space type code.";
479 } else {
480 return loc->second;
481 }
482 }
483
484 bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
485 {
486 FunctionSpaceNamesMapType::iterator loc;
487 loc=m_functionSpaceTypeNames.find(functionSpaceType);
488 return (loc!=m_functionSpaceTypeNames.end());
489 }
490
491 void MeshAdapter::setFunctionSpaceTypeNames()
492 {
493 m_functionSpaceTypeNames.insert
494 (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));
495 m_functionSpaceTypeNames.insert
496 (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));
497 m_functionSpaceTypeNames.insert
498 (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
499 m_functionSpaceTypeNames.insert
500 (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
501 m_functionSpaceTypeNames.insert
502 (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
503 m_functionSpaceTypeNames.insert
504 (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
505 m_functionSpaceTypeNames.insert
506 (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
507 m_functionSpaceTypeNames.insert
508 (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
509 m_functionSpaceTypeNames.insert
510 (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
511 m_functionSpaceTypeNames.insert
512 (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
513 m_functionSpaceTypeNames.insert
514 (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
515 m_functionSpaceTypeNames.insert
516 (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
517 m_functionSpaceTypeNames.insert
518 (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
519 }
520
521 int MeshAdapter::getContinuousFunctionCode() const
522 {
523 return Nodes;
524 }
525 int MeshAdapter::getReducedContinuousFunctionCode() const
526 {
527 return ReducedNodes;
528 }
529
530 int MeshAdapter::getFunctionCode() const
531 {
532 return Elements;
533 }
534 int MeshAdapter::getReducedFunctionCode() const
535 {
536 return ReducedElements;
537 }
538
539 int MeshAdapter::getFunctionOnBoundaryCode() const
540 {
541 return FaceElements;
542 }
543 int MeshAdapter::getReducedFunctionOnBoundaryCode() const
544 {
545 return ReducedFaceElements;
546 }
547
548 int MeshAdapter::getFunctionOnContactZeroCode() const
549 {
550 return ContactElementsZero;
551 }
552 int MeshAdapter::getReducedFunctionOnContactZeroCode() const
553 {
554 return ReducedContactElementsZero;
555 }
556
557 int MeshAdapter::getFunctionOnContactOneCode() const
558 {
559 return ContactElementsOne;
560 }
561 int MeshAdapter::getReducedFunctionOnContactOneCode() const
562 {
563 return ReducedContactElementsOne;
564 }
565
566 int MeshAdapter::getSolutionCode() const
567 {
568 return DegreesOfFreedom;
569 }
570
571 int MeshAdapter::getReducedSolutionCode() const
572 {
573 return ReducedDegreesOfFreedom;
574 }
575
576 int MeshAdapter::getDiracDeltaFunctionCode() const
577 {
578 return Points;
579 }
580
581 //
582 // return the spatial dimension of the Mesh:
583 //
584 int MeshAdapter::getDim() const
585 {
586 int numDim=Finley_Mesh_getDim(m_finleyMesh.get());
587 checkFinleyError();
588 return numDim;
589 }
590
591 //
592 // return the number of data points per sample and the number of samples
593 // needed to represent data on a parts of the mesh.
594 //
595 pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
596 {
597 int numDataPointsPerSample=0;
598 int numSamples=0;
599 Finley_Mesh* mesh=m_finleyMesh.get();
600 switch (functionSpaceCode) {
601 case(Nodes):
602 numDataPointsPerSample=1;
603 numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
604 break;
605 case(ReducedNodes):
606 numDataPointsPerSample=1;
607 numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
608 break;
609 case(Elements):
610 if (mesh->Elements!=NULL) {
611 numSamples=mesh->Elements->numElements;
612 numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
613 }
614 break;
615 case(ReducedElements):
616 if (mesh->Elements!=NULL) {
617 numSamples=mesh->Elements->numElements;
618 numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;
619 }
620 break;
621 case(FaceElements):
622 if (mesh->FaceElements!=NULL) {
623 numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
624 numSamples=mesh->FaceElements->numElements;
625 }
626 break;
627 case(ReducedFaceElements):
628 if (mesh->FaceElements!=NULL) {
629 numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;
630 numSamples=mesh->FaceElements->numElements;
631 }
632 break;
633 case(Points):
634 if (mesh->Points!=NULL) {
635 numDataPointsPerSample=1;
636 numSamples=mesh->Points->numElements;
637 }
638 break;
639 case(ContactElementsZero):
640 if (mesh->ContactElements!=NULL) {
641 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
642 numSamples=mesh->ContactElements->numElements;
643 }
644 break;
645 case(ReducedContactElementsZero):
646 if (mesh->ContactElements!=NULL) {
647 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
648 numSamples=mesh->ContactElements->numElements;
649 }
650 break;
651 case(ContactElementsOne):
652 if (mesh->ContactElements!=NULL) {
653 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
654 numSamples=mesh->ContactElements->numElements;
655 }
656 break;
657 case(ReducedContactElementsOne):
658 if (mesh->ContactElements!=NULL) {
659 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
660 numSamples=mesh->ContactElements->numElements;
661 }
662 break;
663 case(DegreesOfFreedom):
664 if (mesh->Nodes!=NULL) {
665 numDataPointsPerSample=1;
666 numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
667 }
668 break;
669 case(ReducedDegreesOfFreedom):
670 if (mesh->Nodes!=NULL) {
671 numDataPointsPerSample=1;
672 numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
673 }
674 break;
675 default:
676 stringstream temp;
677 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
678 throw FinleyAdapterException(temp.str());
679 break;
680 }
681 return pair<int,int>(numDataPointsPerSample,numSamples);
682 }
683
684 //
685 // adds linear PDE of second order into a given stiffness matrix and right hand side:
686 //
687 void MeshAdapter::addPDEToSystem(
688 SystemMatrixAdapter& mat, escript::Data& rhs,
689 const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,const escript::Data& X,const escript::Data& Y,
690 const escript::Data& d, const escript::Data& y,
691 const escript::Data& d_contact,const escript::Data& y_contact) const
692 {
693 escriptDataC _rhs=rhs.getDataC();
694 escriptDataC _A =A.getDataC();
695 escriptDataC _B=B.getDataC();
696 escriptDataC _C=C.getDataC();
697 escriptDataC _D=D.getDataC();
698 escriptDataC _X=X.getDataC();
699 escriptDataC _Y=Y.getDataC();
700 escriptDataC _d=d.getDataC();
701 escriptDataC _y=y.getDataC();
702 escriptDataC _d_contact=d_contact.getDataC();
703 escriptDataC _y_contact=y_contact.getDataC();
704
705 Finley_Mesh* mesh=m_finleyMesh.get();
706
707 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
708 checkFinleyError();
709
710 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
711 checkFinleyError();
712
713 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
714 checkFinleyError();
715 }
716
717 void MeshAdapter::addPDEToLumpedSystem(
718 escript::Data& mat,
719 const escript::Data& D,
720 const escript::Data& d) const
721 {
722 escriptDataC _mat=mat.getDataC();
723 escriptDataC _D=D.getDataC();
724 escriptDataC _d=d.getDataC();
725
726 Finley_Mesh* mesh=m_finleyMesh.get();
727
728 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
729 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
730
731 checkFinleyError();
732 }
733
734
735 //
736 // adds linear PDE of second order into the right hand side only
737 //
738 void MeshAdapter::addPDEToRHS( escript::Data& rhs, const escript::Data& X,const escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
739 {
740 Finley_Mesh* mesh=m_finleyMesh.get();
741
742 escriptDataC _rhs=rhs.getDataC();
743 escriptDataC _X=X.getDataC();
744 escriptDataC _Y=Y.getDataC();
745 escriptDataC _y=y.getDataC();
746 escriptDataC _y_contact=y_contact.getDataC();
747
748 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
749 checkFinleyError();
750
751 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
752 checkFinleyError();
753
754 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
755 checkFinleyError();
756 }
757
758 //
759 // interpolates data between different function spaces:
760 //
761 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
762 {
763 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
764 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
765 if (inDomain!=*this)
766 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
767 if (targetDomain!=*this)
768 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
769
770 Finley_Mesh* mesh=m_finleyMesh.get();
771 escriptDataC _target=target.getDataC();
772 escriptDataC _in=in.getDataC();
773 switch(in.getFunctionSpace().getTypeCode()) {
774 case(Nodes):
775 switch(target.getFunctionSpace().getTypeCode()) {
776 case(Nodes):
777 case(ReducedNodes):
778 case(DegreesOfFreedom):
779 case(ReducedDegreesOfFreedom):
780 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
781 break;
782 case(Elements):
783 case(ReducedElements):
784 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
785 break;
786 case(FaceElements):
787 case(ReducedFaceElements):
788 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
789 break;
790 case(Points):
791 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
792 break;
793 case(ContactElementsZero):
794 case(ReducedContactElementsZero):
795 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
796 break;
797 case(ContactElementsOne):
798 case(ReducedContactElementsOne):
799 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
800 break;
801 default:
802 stringstream temp;
803 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
804 throw FinleyAdapterException(temp.str());
805 break;
806 }
807 break;
808 case(ReducedNodes):
809 switch(target.getFunctionSpace().getTypeCode()) {
810 case(Nodes):
811 case(ReducedNodes):
812 case(DegreesOfFreedom):
813 case(ReducedDegreesOfFreedom):
814 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
815 break;
816 case(Elements):
817 case(ReducedElements):
818 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
819 break;
820 case(FaceElements):
821 case(ReducedFaceElements):
822 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
823 break;
824 case(Points):
825 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
826 break;
827 case(ContactElementsZero):
828 case(ReducedContactElementsZero):
829 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
830 break;
831 case(ContactElementsOne):
832 case(ReducedContactElementsOne):
833 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
834 break;
835 default:
836 stringstream temp;
837 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
838 throw FinleyAdapterException(temp.str());
839 break;
840 }
841 break;
842 case(Elements):
843 if (target.getFunctionSpace().getTypeCode()==Elements) {
844 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
845 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
846 Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
847 } else {
848 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
849 }
850 break;
851 case(ReducedElements):
852 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
853 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
854 } else {
855 throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
856 }
857 break;
858 case(FaceElements):
859 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
860 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
861 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
862 Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
863 } else {
864 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
865 }
866 break;
867 case(ReducedFaceElements):
868 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
869 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
870 } else {
871 throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
872 }
873 break;
874 case(Points):
875 if (target.getFunctionSpace().getTypeCode()==Points) {
876 Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
877 } else {
878 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
879 }
880 break;
881 case(ContactElementsZero):
882 case(ContactElementsOne):
883 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
884 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
885 } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
886 Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
887 } else {
888 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
889 break;
890 }
891 break;
892 case(ReducedContactElementsZero):
893 case(ReducedContactElementsOne):
894 if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
895 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
896 } else {
897 throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
898 break;
899 }
900 break;
901 case(DegreesOfFreedom):
902 switch(target.getFunctionSpace().getTypeCode()) {
903 case(ReducedDegreesOfFreedom):
904 case(DegreesOfFreedom):
905 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
906 break;
907
908 case(Nodes):
909 case(ReducedNodes):
910 if (getMPISize()>1) {
911 escript::Data temp=escript::Data(in);
912 temp.expand();
913 escriptDataC _in2 = temp.getDataC();
914 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
915 } else {
916 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
917 }
918 break;
919 case(Elements):
920 case(ReducedElements):
921 if (getMPISize()>1) {
922 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
923 escriptDataC _in2 = temp.getDataC();
924 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
925 } else {
926 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
927 }
928 break;
929 case(FaceElements):
930 case(ReducedFaceElements):
931 if (getMPISize()>1) {
932 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
933 escriptDataC _in2 = temp.getDataC();
934 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
935
936 } else {
937 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
938 }
939 break;
940 case(Points):
941 if (getMPISize()>1) {
942 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
943 escriptDataC _in2 = temp.getDataC();
944 } else {
945 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
946 }
947 break;
948 case(ContactElementsZero):
949 case(ContactElementsOne):
950 case(ReducedContactElementsZero):
951 case(ReducedContactElementsOne):
952 if (getMPISize()>1) {
953 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
954 escriptDataC _in2 = temp.getDataC();
955 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
956 } else {
957 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
958 }
959 break;
960 default:
961 stringstream temp;
962 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
963 throw FinleyAdapterException(temp.str());
964 break;
965 }
966 break;
967 case(ReducedDegreesOfFreedom):
968 switch(target.getFunctionSpace().getTypeCode()) {
969 case(Nodes):
970 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
971 break;
972 case(ReducedNodes):
973 if (getMPISize()>1) {
974 escript::Data temp=escript::Data(in);
975 temp.expand();
976 escriptDataC _in2 = temp.getDataC();
977 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
978 } else {
979 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
980 }
981 break;
982 case(DegreesOfFreedom):
983 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
984 break;
985 case(ReducedDegreesOfFreedom):
986 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
987 break;
988 case(Elements):
989 case(ReducedElements):
990 if (getMPISize()>1) {
991 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
992 escriptDataC _in2 = temp.getDataC();
993 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
994 } else {
995 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
996 }
997 break;
998 case(FaceElements):
999 case(ReducedFaceElements):
1000 if (getMPISize()>1) {
1001 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
1002 escriptDataC _in2 = temp.getDataC();
1003 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1004 } else {
1005 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1006 }
1007 break;
1008 case(Points):
1009 if (getMPISize()>1) {
1010 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
1011 escriptDataC _in2 = temp.getDataC();
1012 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1013 } else {
1014 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1015 }
1016 break;
1017 case(ContactElementsZero):
1018 case(ContactElementsOne):
1019 case(ReducedContactElementsZero):
1020 case(ReducedContactElementsOne):
1021 if (getMPISize()>1) {
1022 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
1023 escriptDataC _in2 = temp.getDataC();
1024 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1025 } else {
1026 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1027 }
1028 break;
1029 default:
1030 stringstream temp;
1031 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1032 throw FinleyAdapterException(temp.str());
1033 break;
1034 }
1035 break;
1036 default:
1037 stringstream temp;
1038 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1039 throw FinleyAdapterException(temp.str());
1040 break;
1041 }
1042 checkFinleyError();
1043 }
1044
1045 //
1046 // copies the locations of sample points into x:
1047 //
1048 void MeshAdapter::setToX(escript::Data& arg) const
1049 {
1050 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
1051 if (argDomain!=*this)
1052 throw FinleyAdapterException("Error - Illegal domain of data point locations");
1053 Finley_Mesh* mesh=m_finleyMesh.get();
1054 // in case of values node coordinates we can do the job directly:
1055 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1056 escriptDataC _arg=arg.getDataC();
1057 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1058 } else {
1059 escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
1060 escriptDataC _tmp_data=tmp_data.getDataC();
1061 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1062 // this is then interpolated onto arg:
1063 interpolateOnDomain(arg,tmp_data);
1064 }
1065 checkFinleyError();
1066 }
1067
1068 //
1069 // return the normal vectors at the location of data points as a Data object:
1070 //
1071 void MeshAdapter::setToNormal(escript::Data& normal) const
1072 {
1073 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
1074 if (normalDomain!=*this)
1075 throw FinleyAdapterException("Error - Illegal domain of normal locations");
1076 Finley_Mesh* mesh=m_finleyMesh.get();
1077 escriptDataC _normal=normal.getDataC();
1078 switch(normal.getFunctionSpace().getTypeCode()) {
1079 case(Nodes):
1080 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1081 break;
1082 case(ReducedNodes):
1083 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1084 break;
1085 case(Elements):
1086 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1087 break;
1088 case(ReducedElements):
1089 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1090 break;
1091 case (FaceElements):
1092 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1093 break;
1094 case (ReducedFaceElements):
1095 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1096 break;
1097 case(Points):
1098 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1099 break;
1100 case (ContactElementsOne):
1101 case (ContactElementsZero):
1102 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1103 break;
1104 case (ReducedContactElementsOne):
1105 case (ReducedContactElementsZero):
1106 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1107 break;
1108 case(DegreesOfFreedom):
1109 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1110 break;
1111 case(ReducedDegreesOfFreedom):
1112 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1113 break;
1114 default:
1115 stringstream temp;
1116 temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1117 throw FinleyAdapterException(temp.str());
1118 break;
1119 }
1120 checkFinleyError();
1121 }
1122
1123 //
1124 // interpolates data to other domain:
1125 //
1126 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1127 {
1128 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
1129 if (targetDomain!=*this)
1130 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1131
1132 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
1133 }
1134
1135 //
1136 // calculates the integral of a function defined of arg:
1137 //
1138 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
1139 {
1140 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
1141 if (argDomain!=*this)
1142 throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1143
1144 double blocktimer_start = blocktimer_time();
1145 Finley_Mesh* mesh=m_finleyMesh.get();
1146 escriptDataC _temp;
1147 escript::Data temp;
1148 escriptDataC _arg=arg.getDataC();
1149 switch(arg.getFunctionSpace().getTypeCode()) {
1150 case(Nodes):
1151 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
1152 _temp=temp.getDataC();
1153 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1154 break;
1155 case(ReducedNodes):
1156 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
1157 _temp=temp.getDataC();
1158 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1159 break;
1160 case(Elements):
1161 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1162 break;
1163 case(ReducedElements):
1164 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1165 break;
1166 case(FaceElements):
1167 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1168 break;
1169 case(ReducedFaceElements):
1170 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1171 break;
1172 case(Points):
1173 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1174 break;
1175 case(ContactElementsZero):
1176 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1177 break;
1178 case(ReducedContactElementsZero):
1179 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1180 break;
1181 case(ContactElementsOne):
1182 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1183 break;
1184 case(ReducedContactElementsOne):
1185 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1186 break;
1187 case(DegreesOfFreedom):
1188 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
1189 _temp=temp.getDataC();
1190 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1191 break;
1192 case(ReducedDegreesOfFreedom):
1193 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
1194 _temp=temp.getDataC();
1195 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1196 break;
1197 default:
1198 stringstream temp;
1199 temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1200 throw FinleyAdapterException(temp.str());
1201 break;
1202 }
1203 checkFinleyError();
1204 blocktimer_increment("integrate()", blocktimer_start);
1205 }
1206
1207 //
1208 // calculates the gradient of arg:
1209 //
1210 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1211 {
1212 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
1213 if (argDomain!=*this)
1214 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1215 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
1216 if (gradDomain!=*this)
1217 throw FinleyAdapterException("Error - Illegal domain of gradient");
1218
1219 Finley_Mesh* mesh=m_finleyMesh.get();
1220 escriptDataC _grad=grad.getDataC();
1221 escriptDataC nodeDataC;
1222 escript::Data temp;
1223 if (getMPISize()>1) {
1224 if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1225 temp=escript::Data( arg, continuousFunction(asAbstractContinuousDomain()) );
1226 nodeDataC = temp.getDataC();
1227 } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1228 temp=escript::Data( arg, reducedContinuousFunction(asAbstractContinuousDomain()) );
1229 nodeDataC = temp.getDataC();
1230 } else {
1231 nodeDataC = arg.getDataC();
1232 }
1233 } else {
1234 nodeDataC = arg.getDataC();
1235 }
1236 switch(grad.getFunctionSpace().getTypeCode()) {
1237 case(Nodes):
1238 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1239 break;
1240 case(ReducedNodes):
1241 throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1242 break;
1243 case(Elements):
1244 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1245 break;
1246 case(ReducedElements):
1247 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1248 break;
1249 case(FaceElements):
1250 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1251 break;
1252 case(ReducedFaceElements):
1253 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1254 break;
1255 case(Points):
1256 throw FinleyAdapterException("Error - Gradient at points is not supported.");
1257 break;
1258 case(ContactElementsZero):
1259 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1260 break;
1261 case(ReducedContactElementsZero):
1262 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1263 break;
1264 case(ContactElementsOne):
1265 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1266 break;
1267 case(ReducedContactElementsOne):
1268 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1269 break;
1270 case(DegreesOfFreedom):
1271 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1272 break;
1273 case(ReducedDegreesOfFreedom):
1274 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1275 break;
1276 default:
1277 stringstream temp;
1278 temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1279 throw FinleyAdapterException(temp.str());
1280 break;
1281 }
1282 checkFinleyError();
1283 }
1284
1285 //
1286 // returns the size of elements:
1287 //
1288 void MeshAdapter::setToSize(escript::Data& size) const
1289 {
1290 Finley_Mesh* mesh=m_finleyMesh.get();
1291 escriptDataC tmp=size.getDataC();
1292 switch(size.getFunctionSpace().getTypeCode()) {
1293 case(Nodes):
1294 throw FinleyAdapterException("Error - Size of nodes is not supported.");
1295 break;
1296 case(ReducedNodes):
1297 throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1298 break;
1299 case(Elements):
1300 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1301 break;
1302 case(ReducedElements):
1303 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1304 break;
1305 case(FaceElements):
1306 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1307 break;
1308 case(ReducedFaceElements):
1309 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1310 break;
1311 case(Points):
1312 throw FinleyAdapterException("Error - Size of point elements is not supported.");
1313 break;
1314 case(ContactElementsZero):
1315 case(ContactElementsOne):
1316 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1317 break;
1318 case(ReducedContactElementsZero):
1319 case(ReducedContactElementsOne):
1320 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1321 break;
1322 case(DegreesOfFreedom):
1323 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1324 break;
1325 case(ReducedDegreesOfFreedom):
1326 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1327 break;
1328 default:
1329 stringstream temp;
1330 temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1331 throw FinleyAdapterException(temp.str());
1332 break;
1333 }
1334 checkFinleyError();
1335 }
1336
1337 // sets the location of nodes:
1338 void MeshAdapter::setNewX(const escript::Data& new_x)
1339 {
1340 Finley_Mesh* mesh=m_finleyMesh.get();
1341 escriptDataC tmp;
1342 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
1343 if (newDomain!=*this)
1344 throw FinleyAdapterException("Error - Illegal domain of new point locations");
1345 tmp = new_x.getDataC();
1346 Finley_Mesh_setCoordinates(mesh,&tmp);
1347 checkFinleyError();
1348 }
1349
1350 // saves a data array in openDX format:
1351 void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
1352 {
1353 int MAX_namelength=256;
1354 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1355 /* win32 refactor */
1356 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1357 for(int i=0;i<num_data;i++)
1358 {
1359 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
1360 }
1361
1362 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1363 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
1364 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
1365
1366 boost::python::list keys=arg.keys();
1367 for (int i=0;i<num_data;++i) {
1368 std::string n=boost::python::extract<std::string>(keys[i]);
1369 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1370 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1371 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
1372 data[i]=d.getDataC();
1373 ptr_data[i]=&(data[i]);
1374 c_names[i]=&(names[i][0]);
1375 if (n.length()>MAX_namelength-1) {
1376 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
1377 c_names[i][MAX_namelength-1]='\0';
1378 } else {
1379 strcpy(c_names[i],n.c_str());
1380 }
1381 }
1382 Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1383 checkFinleyError();
1384
1385 /* win32 refactor */
1386 TMPMEMFREE(c_names);
1387 TMPMEMFREE(data);
1388 TMPMEMFREE(ptr_data);
1389 for(int i=0;i<num_data;i++)
1390 {
1391 TMPMEMFREE(names[i]);
1392 }
1393 TMPMEMFREE(names);
1394
1395 return;
1396 }
1397
1398 // saves a data array in openVTK format:
1399 void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
1400 {
1401 int MAX_namelength=256;
1402 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1403 /* win32 refactor */
1404 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1405 for(int i=0;i<num_data;i++)
1406 {
1407 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
1408 }
1409
1410 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1411 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
1412 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
1413
1414 boost::python::list keys=arg.keys();
1415 for (int i=0;i<num_data;++i) {
1416 std::string n=boost::python::extract<std::string>(keys[i]);
1417 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1418 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1419 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
1420 data[i]=d.getDataC();
1421 ptr_data[i]=&(data[i]);
1422 c_names[i]=&(names[i][0]);
1423 if (n.length()>MAX_namelength-1) {
1424 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
1425 c_names[i][MAX_namelength-1]='\0';
1426 } else {
1427 strcpy(c_names[i],n.c_str());
1428 }
1429 }
1430 Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1431
1432 checkFinleyError();
1433 /* win32 refactor */
1434 TMPMEMFREE(c_names);
1435 TMPMEMFREE(data);
1436 TMPMEMFREE(ptr_data);
1437 for(int i=0;i<num_data;i++)
1438 {
1439 TMPMEMFREE(names[i]);
1440 }
1441 TMPMEMFREE(names);
1442
1443 return;
1444 }
1445
1446
1447 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
1448 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1449 const int row_blocksize,
1450 const escript::FunctionSpace& row_functionspace,
1451 const int column_blocksize,
1452 const escript::FunctionSpace& column_functionspace,
1453 const int type) const
1454 {
1455 int reduceRowOrder=0;
1456 int reduceColOrder=0;
1457 // is the domain right?
1458 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
1459 if (row_domain!=*this)
1460 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1461 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
1462 if (col_domain!=*this)
1463 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1464 // is the function space type right
1465 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1466 reduceRowOrder=0;
1467 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1468 reduceRowOrder=1;
1469 } else {
1470 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1471 }
1472 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1473 reduceColOrder=0;
1474 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1475 reduceColOrder=1;
1476 } else {
1477 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1478 }
1479 // generate matrix:
1480
1481 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1482 checkFinleyError();
1483 Paso_SystemMatrix* fsystemMatrix;
1484 int trilinos = 0;
1485 if (trilinos) {
1486 #ifdef TRILINOS
1487 /* Allocation Epetra_VrbMatrix here */
1488 #endif
1489 }
1490 else {
1491 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
1492 }
1493 checkPasoError();
1494 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1495 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1496 }
1497
1498 //
1499 // vtkObject MeshAdapter::createVtkObject() const
1500 // TODO:
1501 //
1502
1503 //
1504 // returns true if data at the atom_type is considered as being cell centered:
1505 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1506 {
1507 switch(functionSpaceCode) {
1508 case(Nodes):
1509 case(DegreesOfFreedom):
1510 case(ReducedDegreesOfFreedom):
1511 return false;
1512 break;
1513 case(Elements):
1514 case(FaceElements):
1515 case(Points):
1516 case(ContactElementsZero):
1517 case(ContactElementsOne):
1518 case(ReducedElements):
1519 case(ReducedFaceElements):
1520 case(ReducedContactElementsZero):
1521 case(ReducedContactElementsOne):
1522 return true;
1523 break;
1524 default:
1525 stringstream temp;
1526 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1527 throw FinleyAdapterException(temp.str());
1528 break;
1529 }
1530 checkFinleyError();
1531 return false;
1532 }
1533
1534 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1535 {
1536 switch(functionSpaceType_source) {
1537 case(Nodes):
1538 switch(functionSpaceType_target) {
1539 case(Nodes):
1540 case(ReducedNodes):
1541 case(ReducedDegreesOfFreedom):
1542 case(DegreesOfFreedom):
1543 case(Elements):
1544 case(ReducedElements):
1545 case(FaceElements):
1546 case(ReducedFaceElements):
1547 case(Points):
1548 case(ContactElementsZero):
1549 case(ReducedContactElementsZero):
1550 case(ContactElementsOne):
1551 case(ReducedContactElementsOne):
1552 return true;
1553 default:
1554 stringstream temp;
1555 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1556 throw FinleyAdapterException(temp.str());
1557 }
1558 break;
1559 case(ReducedNodes):
1560 switch(functionSpaceType_target) {
1561 case(ReducedNodes):
1562 case(ReducedDegreesOfFreedom):
1563 case(Elements):
1564 case(ReducedElements):
1565 case(FaceElements):
1566 case(ReducedFaceElements):
1567 case(Points):
1568 case(ContactElementsZero):
1569 case(ReducedContactElementsZero):
1570 case(ContactElementsOne):
1571 case(ReducedContactElementsOne):
1572 return true;
1573 case(Nodes):
1574 case(DegreesOfFreedom):
1575 return false;
1576 default:
1577 stringstream temp;
1578 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1579 throw FinleyAdapterException(temp.str());
1580 }
1581 break;
1582 case(Elements):
1583 if (functionSpaceType_target==Elements) {
1584 return true;
1585 } else if (functionSpaceType_target==ReducedElements) {
1586 return true;
1587 } else {
1588 return false;
1589 }
1590 case(ReducedElements):
1591 if (functionSpaceType_target==ReducedElements) {
1592 return true;
1593 } else {
1594 return false;
1595 }
1596 case(FaceElements):
1597 if (functionSpaceType_target==FaceElements) {
1598 return true;
1599 } else if (functionSpaceType_target==ReducedFaceElements) {
1600 return true;
1601 } else {
1602 return false;
1603 }
1604 case(ReducedFaceElements):
1605 if (functionSpaceType_target==ReducedFaceElements) {
1606 return true;
1607 } else {
1608 return false;
1609 }
1610 case(Points):
1611 if (functionSpaceType_target==Points) {
1612 return true;
1613 } else {
1614 return false;
1615 }
1616 case(ContactElementsZero):
1617 case(ContactElementsOne):
1618 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1619 return true;
1620 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1621 return true;
1622 } else {
1623 return false;
1624 }
1625 case(ReducedContactElementsZero):
1626 case(ReducedContactElementsOne):
1627 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1628 return true;
1629 } else {
1630 return false;
1631 }
1632 case(DegreesOfFreedom):
1633 switch(functionSpaceType_target) {
1634 case(ReducedDegreesOfFreedom):
1635 case(DegreesOfFreedom):
1636 case(Nodes):
1637 case(ReducedNodes):
1638 case(Elements):
1639 case(ReducedElements):
1640 case(Points):
1641 case(FaceElements):
1642 case(ReducedFaceElements):
1643 case(ContactElementsZero):
1644 case(ReducedContactElementsZero):
1645 case(ContactElementsOne):
1646 case(ReducedContactElementsOne):
1647 return true;
1648 default:
1649 stringstream temp;
1650 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1651 throw FinleyAdapterException(temp.str());
1652 }
1653 break;
1654 case(ReducedDegreesOfFreedom):
1655 switch(functionSpaceType_target) {
1656 case(ReducedDegreesOfFreedom):
1657 case(ReducedNodes):
1658 case(Elements):
1659 case(ReducedElements):
1660 case(FaceElements):
1661 case(ReducedFaceElements):
1662 case(Points):
1663 case(ContactElementsZero):
1664 case(ReducedContactElementsZero):
1665 case(ContactElementsOne):
1666 case(ReducedContactElementsOne):
1667 return true;
1668 case(Nodes):
1669 case(DegreesOfFreedom):
1670 return false;
1671 default:
1672 stringstream temp;
1673 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1674 throw FinleyAdapterException(temp.str());
1675 }
1676 break;
1677 default:
1678 stringstream temp;
1679 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1680 throw FinleyAdapterException(temp.str());
1681 break;
1682 }
1683 checkFinleyError();
1684 return false;
1685 }
1686
1687 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1688 {
1689 return false;
1690 }
1691
1692 bool MeshAdapter::operator==(const AbstractDomain& other) const
1693 {
1694 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1695 if (temp!=0) {
1696 return (m_finleyMesh==temp->m_finleyMesh);
1697 } else {
1698 return false;
1699 }
1700 }
1701
1702 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1703 {
1704 return !(operator==(other));
1705 }
1706
1707 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1708 {
1709 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1710 checkPasoError();
1711 return out;
1712 }
1713
1714 escript::Data MeshAdapter::getX() const
1715 {
1716 return continuousFunction(asAbstractContinuousDomain()).getX();
1717 }
1718
1719 escript::Data MeshAdapter::getNormal() const
1720 {
1721 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1722 }
1723
1724 escript::Data MeshAdapter::getSize() const
1725 {
1726 return function(asAbstractContinuousDomain()).getSize();
1727 }
1728
1729 int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1730 {
1731 int *out=0,i;
1732 Finley_Mesh* mesh=m_finleyMesh.get();
1733 switch (functionSpaceType) {
1734 case(Nodes):
1735 out=mesh->Nodes->Id;
1736 break;
1737 case(ReducedNodes):
1738 out=mesh->Nodes->reducedNodesId;
1739 break;
1740 case(Elements):
1741 out=mesh->Elements->Id;
1742 break;
1743 case(ReducedElements):
1744 out=mesh->Elements->Id;
1745 break;
1746 case(FaceElements):
1747 out=mesh->FaceElements->Id;
1748 break;
1749 case(ReducedFaceElements):
1750 out=mesh->FaceElements->Id;
1751 break;
1752 case(Points):
1753 out=mesh->Points->Id;
1754 break;
1755 case(ContactElementsZero):
1756 out=mesh->ContactElements->Id;
1757 break;
1758 case(ReducedContactElementsZero):
1759 out=mesh->ContactElements->Id;
1760 break;
1761 case(ContactElementsOne):
1762 out=mesh->ContactElements->Id;
1763 break;
1764 case(ReducedContactElementsOne):
1765 out=mesh->ContactElements->Id;
1766 break;
1767 case(DegreesOfFreedom):
1768 out=mesh->Nodes->degreesOfFreedomId;
1769 break;
1770 case(ReducedDegreesOfFreedom):
1771 out=mesh->Nodes->reducedDegreesOfFreedomId;
1772 break;
1773 default:
1774 stringstream temp;
1775 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1776 throw FinleyAdapterException(temp.str());
1777 break;
1778 }
1779 return out;
1780 }
1781 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1782 {
1783 int out=0;
1784 Finley_Mesh* mesh=m_finleyMesh.get();
1785 switch (functionSpaceType) {
1786 case(Nodes):
1787 out=mesh->Nodes->Tag[sampleNo];
1788 break;
1789 case(ReducedNodes):
1790 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1791 break;
1792 case(Elements):
1793 out=mesh->Elements->Tag[sampleNo];
1794 break;
1795 case(ReducedElements):
1796 out=mesh->Elements->Tag[sampleNo];
1797 break;
1798 case(FaceElements):
1799 out=mesh->FaceElements->Tag[sampleNo];
1800 break;
1801 case(ReducedFaceElements):
1802 out=mesh->FaceElements->Tag[sampleNo];
1803 break;
1804 case(Points):
1805 out=mesh->Points->Tag[sampleNo];
1806 break;
1807 case(ContactElementsZero):
1808 out=mesh->ContactElements->Tag[sampleNo];
1809 break;
1810 case(ReducedContactElementsZero):
1811 out=mesh->ContactElements->Tag[sampleNo];
1812 break;
1813 case(ContactElementsOne):
1814 out=mesh->ContactElements->Tag[sampleNo];
1815 break;
1816 case(ReducedContactElementsOne):
1817 out=mesh->ContactElements->Tag[sampleNo];
1818 break;
1819 case(DegreesOfFreedom):
1820 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1821 break;
1822 case(ReducedDegreesOfFreedom):
1823 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1824 break;
1825 default:
1826 stringstream temp;
1827 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1828 throw FinleyAdapterException(temp.str());
1829 break;
1830 }
1831 return out;
1832 }
1833
1834
1835 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1836 {
1837 Finley_Mesh* mesh=m_finleyMesh.get();
1838 escriptDataC tmp=mask.getDataC();
1839 switch(functionSpaceType) {
1840 case(Nodes):
1841 Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1842 break;
1843 case(ReducedNodes):
1844 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
1845 break;
1846 case(DegreesOfFreedom):
1847 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1848 break;
1849 case(ReducedDegreesOfFreedom):
1850 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1851 break;
1852 case(Elements):
1853 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1854 break;
1855 case(ReducedElements):
1856 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1857 break;
1858 case(FaceElements):
1859 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1860 break;
1861 case(ReducedFaceElements):
1862 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1863 break;
1864 case(Points):
1865 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1866 break;
1867 case(ContactElementsZero):
1868 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1869 break;
1870 case(ReducedContactElementsZero):
1871 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1872 break;
1873 case(ContactElementsOne):
1874 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1875 break;
1876 case(ReducedContactElementsOne):
1877 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1878 break;
1879 default:
1880 stringstream temp;
1881 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1882 throw FinleyAdapterException(temp.str());
1883 }
1884 checkFinleyError();
1885 return;
1886 }
1887
1888 void MeshAdapter::setTagMap(const std::string& name, int tag)
1889 {
1890 Finley_Mesh* mesh=m_finleyMesh.get();
1891 Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
1892 checkPasoError();
1893 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1894 }
1895
1896 int MeshAdapter::getTag(const std::string& name) const
1897 {
1898 Finley_Mesh* mesh=m_finleyMesh.get();
1899 int tag=0;
1900 tag=Finley_Mesh_getTag(mesh, name.c_str());
1901 checkPasoError();
1902 // throwStandardException("MeshAdapter::getTag is not implemented.");
1903 return tag;
1904 }
1905
1906 bool MeshAdapter::isValidTagName(const std::string& name) const
1907 {
1908 Finley_Mesh* mesh=m_finleyMesh.get();
1909 return Finley_Mesh_isValidTagName(mesh,name.c_str());
1910 }
1911
1912 std::string MeshAdapter::showTagNames() const
1913 {
1914 stringstream temp;
1915 Finley_Mesh* mesh=m_finleyMesh.get();
1916 Finley_TagMap* tag_map=mesh->TagMap;
1917 while (tag_map) {
1918 temp << tag_map->name;
1919 tag_map=tag_map->next;
1920 if (tag_map) temp << ", ";
1921 }
1922 return temp.str();
1923 }
1924
1925 } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26