/[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 1344 - (show annotations)
Wed Nov 14 04:28:25 2007 UTC (11 years, 10 months ago) by ksteube
File size: 79785 byte(s)
Now using TMPMEMFREE but still thrashes memory for very large mesh.
Future solution will be to use NetCDF subsetting and less temp memory.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26