/[escript]/branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Contents of /branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26