/[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 3231 - (show annotations)
Fri Oct 1 01:53:46 2010 UTC (8 years, 6 months ago) by jfenwick
File size: 72127 byte(s)
Fix MPI and OMP problems not detected in serial

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26