/[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 3090 - (show annotations)
Wed Aug 11 00:51:28 2010 UTC (8 years, 8 months ago) by jfenwick
File size: 72859 byte(s)
Removed:
DUDLEY_CONTACT_ELEMENTS_1
DUDLEY_CONTACT_ELEMENTS_2
DUDLEY_REDUCED_CONTACT_ELEMENTS_1
DUDLEY_REDUCED_CONTACT_ELEMENTS_2

escript tests now query if the domain supports contact elements
before trying to use them.


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->Parametrization->numQuadNodes;
634 }
635 break;
636 case(ReducedElements):
637 if (mesh->Elements!=NULL) {
638 numSamples=mesh->Elements->numElements;
639 numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
640 }
641 break;
642 case(FaceElements):
643 if (mesh->FaceElements!=NULL) {
644 numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
645 numSamples=mesh->FaceElements->numElements;
646 }
647 break;
648 case(ReducedFaceElements):
649 if (mesh->FaceElements!=NULL) {
650 numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->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,
688 const escript::Data& d_contact,const escript::Data& y_contact) const
689 {
690 escriptDataC _rhs=rhs.getDataC();
691 escriptDataC _A =A.getDataC();
692 escriptDataC _B=B.getDataC();
693 escriptDataC _C=C.getDataC();
694 escriptDataC _D=D.getDataC();
695 escriptDataC _X=X.getDataC();
696 escriptDataC _Y=Y.getDataC();
697 escriptDataC _d=d.getDataC();
698 escriptDataC _y=y.getDataC();
699 escriptDataC _d_contact=d_contact.getDataC();
700 escriptDataC _y_contact=y_contact.getDataC();
701
702 Dudley_Mesh* mesh=m_dudleyMesh.get();
703
704 Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
705 checkDudleyError();
706
707 Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
708 checkDudleyError();
709
710 checkDudleyError();
711 }
712
713 void MeshAdapter::addPDEToLumpedSystem(
714 escript::Data& mat,
715 const escript::Data& D,
716 const escript::Data& d) const
717 {
718 escriptDataC _mat=mat.getDataC();
719 escriptDataC _D=D.getDataC();
720 escriptDataC _d=d.getDataC();
721
722 Dudley_Mesh* mesh=m_dudleyMesh.get();
723
724 Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
725 checkDudleyError();
726
727 Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
728 checkDudleyError();
729
730 }
731
732
733 //
734 // adds linear PDE of second order into the right hand side only
735 //
736 void MeshAdapter::addPDEToRHS( escript::Data& rhs, const escript::Data& X,const escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
737 {
738 Dudley_Mesh* mesh=m_dudleyMesh.get();
739
740 escriptDataC _rhs=rhs.getDataC();
741 escriptDataC _X=X.getDataC();
742 escriptDataC _Y=Y.getDataC();
743 escriptDataC _y=y.getDataC();
744 escriptDataC _y_contact=y_contact.getDataC();
745
746 Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
747 checkDudleyError();
748
749 Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
750 checkDudleyError();
751
752 checkDudleyError();
753 }
754 //
755 // adds PDE of second order into a transport problem
756 //
757 void MeshAdapter::addPDEToTransportProblem(
758 TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,
759 const escript::Data& A, const escript::Data& B, const escript::Data& C,
760 const escript::Data& D,const escript::Data& X,const escript::Data& Y,
761 const escript::Data& d, const escript::Data& y,
762 const escript::Data& d_contact,const escript::Data& y_contact) const
763 {
764 DataTypes::ShapeType shape;
765 source.expand();
766 escriptDataC _source=source.getDataC();
767 escriptDataC _M=M.getDataC();
768 escriptDataC _A=A.getDataC();
769 escriptDataC _B=B.getDataC();
770 escriptDataC _C=C.getDataC();
771 escriptDataC _D=D.getDataC();
772 escriptDataC _X=X.getDataC();
773 escriptDataC _Y=Y.getDataC();
774 escriptDataC _d=d.getDataC();
775 escriptDataC _y=y.getDataC();
776 escriptDataC _d_contact=d_contact.getDataC();
777 escriptDataC _y_contact=y_contact.getDataC();
778
779 Dudley_Mesh* mesh=m_dudleyMesh.get();
780 Paso_TransportProblem* _tp = tp.getPaso_TransportProblem();
781
782 Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
783 checkDudleyError();
784
785 Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
786 checkDudleyError();
787
788 Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
789 checkDudleyError();
790
791 checkDudleyError();
792 }
793
794 //
795 // interpolates data between different function spaces:
796 //
797 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
798 {
799 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
800 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
801 if (inDomain!=*this)
802 throw DudleyAdapterException("Error - Illegal domain of interpolant.");
803 if (targetDomain!=*this)
804 throw DudleyAdapterException("Error - Illegal domain of interpolation target.");
805
806 Dudley_Mesh* mesh=m_dudleyMesh.get();
807 escriptDataC _target=target.getDataC();
808 escriptDataC _in=in.getDataC();
809 switch(in.getFunctionSpace().getTypeCode()) {
810 case(Nodes):
811 switch(target.getFunctionSpace().getTypeCode()) {
812 case(Nodes):
813 case(ReducedNodes):
814 case(DegreesOfFreedom):
815 case(ReducedDegreesOfFreedom):
816 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
817 break;
818 case(Elements):
819 case(ReducedElements):
820 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
821 break;
822 case(FaceElements):
823 case(ReducedFaceElements):
824 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
825 break;
826 case(Points):
827 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
828 break;
829 default:
830 stringstream temp;
831 temp << "Error - Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
832 throw DudleyAdapterException(temp.str());
833 break;
834 }
835 break;
836 case(ReducedNodes):
837 switch(target.getFunctionSpace().getTypeCode()) {
838 case(Nodes):
839 case(ReducedNodes):
840 case(DegreesOfFreedom):
841 case(ReducedDegreesOfFreedom):
842 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
843 break;
844 case(Elements):
845 case(ReducedElements):
846 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
847 break;
848 case(FaceElements):
849 case(ReducedFaceElements):
850 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
851 break;
852 case(Points):
853 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
854 break;
855 default:
856 stringstream temp;
857 temp << "Error - Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
858 throw DudleyAdapterException(temp.str());
859 break;
860 }
861 break;
862 case(Elements):
863 if (target.getFunctionSpace().getTypeCode()==Elements) {
864 Dudley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
865 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
866 Dudley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
867 } else {
868 throw DudleyAdapterException("Error - No interpolation with data on elements possible.");
869 }
870 break;
871 case(ReducedElements):
872 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
873 Dudley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
874 } else {
875 throw DudleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
876 }
877 break;
878 case(FaceElements):
879 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
880 Dudley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
881 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
882 Dudley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
883 } else {
884 throw DudleyAdapterException("Error - No interpolation with data on face elements possible.");
885 }
886 break;
887 case(ReducedFaceElements):
888 if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
889 Dudley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
890 } else {
891 throw DudleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
892 }
893 break;
894 case(Points):
895 if (target.getFunctionSpace().getTypeCode()==Points) {
896 Dudley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
897 } else {
898 throw DudleyAdapterException("Error - No interpolation with data on points possible.");
899 }
900 break;
901 case(DegreesOfFreedom):
902 switch(target.getFunctionSpace().getTypeCode()) {
903 case(ReducedDegreesOfFreedom):
904 case(DegreesOfFreedom):
905 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
906 break;
907
908 case(Nodes):
909 case(ReducedNodes):
910 if (getMPISize()>1) {
911 escript::Data temp=escript::Data(in);
912 temp.expand();
913 escriptDataC _in2 = temp.getDataC();
914 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
915 } else {
916 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
917 }
918 break;
919 case(Elements):
920 case(ReducedElements):
921 if (getMPISize()>1) {
922 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
923 escriptDataC _in2 = temp.getDataC();
924 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
925 } else {
926 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
927 }
928 break;
929 case(FaceElements):
930 case(ReducedFaceElements):
931 if (getMPISize()>1) {
932 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
933 escriptDataC _in2 = temp.getDataC();
934 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
935
936 } else {
937 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
938 }
939 break;
940 case(Points):
941 if (getMPISize()>1) {
942 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
943 escriptDataC _in2 = temp.getDataC();
944 } else {
945 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
946 }
947 break;
948 default:
949 stringstream temp;
950 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
951 throw DudleyAdapterException(temp.str());
952 break;
953 }
954 break;
955 case(ReducedDegreesOfFreedom):
956 switch(target.getFunctionSpace().getTypeCode()) {
957 case(Nodes):
958 throw DudleyAdapterException("Error - Dudley does not support interpolation from reduced degrees of freedom to mesh nodes.");
959 break;
960 case(ReducedNodes):
961 if (getMPISize()>1) {
962 escript::Data temp=escript::Data(in);
963 temp.expand();
964 escriptDataC _in2 = temp.getDataC();
965 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
966 } else {
967 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
968 }
969 break;
970 case(DegreesOfFreedom):
971 throw DudleyAdapterException("Error - Dudley does not support interpolation from reduced degrees of freedom to degrees of freedom");
972 break;
973 case(ReducedDegreesOfFreedom):
974 Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
975 break;
976 case(Elements):
977 case(ReducedElements):
978 if (getMPISize()>1) {
979 escript::Data temp=escript::Data( in, reducedContinuousFunction(asAbstractContinuousDomain()) );
980 escriptDataC _in2 = temp.getDataC();
981 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
982 } else {
983 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
984 }
985 break;
986 case(FaceElements):
987 case(ReducedFaceElements):
988 if (getMPISize()>1) {
989 escript::Data temp=escript::Data( in, reducedContinuousFunction(asAbstractContinuousDomain()) );
990 escriptDataC _in2 = temp.getDataC();
991 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
992 } else {
993 Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
994 }
995 break;
996 case(Points):
997 if (getMPISize()>1) {
998 escript::Data temp=escript::Data( in, reducedContinuousFunction(asAbstractContinuousDomain()) );
999 escriptDataC _in2 = temp.getDataC();
1000 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1001 } else {
1002 Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1003 }
1004 break;
1005 default:
1006 stringstream temp;
1007 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1008 throw DudleyAdapterException(temp.str());
1009 break;
1010 }
1011 break;
1012 default:
1013 stringstream temp;
1014 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1015 throw DudleyAdapterException(temp.str());
1016 break;
1017 }
1018 checkDudleyError();
1019 }
1020
1021 //
1022 // copies the locations of sample points into x:
1023 //
1024 void MeshAdapter::setToX(escript::Data& arg) const
1025 {
1026 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1027 if (argDomain!=*this)
1028 throw DudleyAdapterException("Error - Illegal domain of data point locations");
1029 Dudley_Mesh* mesh=m_dudleyMesh.get();
1030 // in case of values node coordinates we can do the job directly:
1031 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1032 escriptDataC _arg=arg.getDataC();
1033 Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1034 } else {
1035 escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
1036 escriptDataC _tmp_data=tmp_data.getDataC();
1037 Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1038 // this is then interpolated onto arg:
1039 interpolateOnDomain(arg,tmp_data);
1040 }
1041 checkDudleyError();
1042 }
1043
1044 //
1045 // return the normal vectors at the location of data points as a Data object:
1046 //
1047 void MeshAdapter::setToNormal(escript::Data& normal) const
1048 {
1049 /* const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1050 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1051 if (normalDomain!=*this)
1052 throw DudleyAdapterException("Error - Illegal domain of normal locations");
1053 Dudley_Mesh* mesh=m_dudleyMesh.get();
1054 escriptDataC _normal=normal.getDataC();
1055 switch(normal.getFunctionSpace().getTypeCode()) {
1056 case(Nodes):
1057 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for nodes");
1058 break;
1059 case(ReducedNodes):
1060 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for reduced nodes");
1061 break;
1062 case(Elements):
1063 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for elements");
1064 break;
1065 case(ReducedElements):
1066 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for elements with reduced integration order");
1067 break;
1068 case (FaceElements):
1069 Dudley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1070 break;
1071 case (ReducedFaceElements):
1072 Dudley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1073 break;
1074 case(Points):
1075 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for point elements");
1076 break;
1077 case(DegreesOfFreedom):
1078 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for degrees of freedom.");
1079 break;
1080 case(ReducedDegreesOfFreedom):
1081 throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for reduced degrees of freedom.");
1082 break;
1083 default:
1084 stringstream temp;
1085 temp << "Error - Normal Vectors: Dudley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1086 throw DudleyAdapterException(temp.str());
1087 break;
1088 }
1089 checkDudleyError();
1090 }
1091
1092 //
1093 // interpolates data to other domain:
1094 //
1095 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1096 {
1097 const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1098 const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1099 if (targetDomain!=this)
1100 throw DudleyAdapterException("Error - Illegal domain of interpolation target");
1101
1102 throw DudleyAdapterException("Error - Dudley does not allow interpolation across domains yet.");
1103 }
1104
1105 //
1106 // calculates the integral of a function defined of arg:
1107 //
1108 void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1109 {
1110 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1111 if (argDomain!=*this)
1112 throw DudleyAdapterException("Error - Illegal domain of integration kernel");
1113
1114 double blocktimer_start = blocktimer_time();
1115 Dudley_Mesh* mesh=m_dudleyMesh.get();
1116 escriptDataC _temp;
1117 escript::Data temp;
1118 escriptDataC _arg=arg.getDataC();
1119 switch(arg.getFunctionSpace().getTypeCode()) {
1120 case(Nodes):
1121 temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1122 _temp=temp.getDataC();
1123 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1124 break;
1125 case(ReducedNodes):
1126 temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1127 _temp=temp.getDataC();
1128 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1129 break;
1130 case(Elements):
1131 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1132 break;
1133 case(ReducedElements):
1134 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1135 break;
1136 case(FaceElements):
1137 Dudley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1138 break;
1139 case(ReducedFaceElements):
1140 Dudley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1141 break;
1142 case(Points):
1143 throw DudleyAdapterException("Error - Integral of data on points is not supported.");
1144 break;
1145 case(DegreesOfFreedom):
1146 temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1147 _temp=temp.getDataC();
1148 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1149 break;
1150 case(ReducedDegreesOfFreedom):
1151 temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1152 _temp=temp.getDataC();
1153 Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1154 break;
1155 default:
1156 stringstream temp;
1157 temp << "Error - Integrals: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1158 throw DudleyAdapterException(temp.str());
1159 break;
1160 }
1161 checkDudleyError();
1162 blocktimer_increment("integrate()", blocktimer_start);
1163 }
1164
1165 //
1166 // calculates the gradient of arg:
1167 //
1168 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1169 {
1170 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1171 if (argDomain!=*this)
1172 throw DudleyAdapterException("Error - Illegal domain of gradient argument");
1173 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1174 if (gradDomain!=*this)
1175 throw DudleyAdapterException("Error - Illegal domain of gradient");
1176
1177 Dudley_Mesh* mesh=m_dudleyMesh.get();
1178 escriptDataC _grad=grad.getDataC();
1179 escriptDataC nodeDataC;
1180 escript::Data temp;
1181 if (getMPISize()>1) {
1182 if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1183 temp=escript::Data( arg, continuousFunction(asAbstractContinuousDomain()) );
1184 nodeDataC = temp.getDataC();
1185 } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1186 temp=escript::Data( arg, reducedContinuousFunction(asAbstractContinuousDomain()) );
1187 nodeDataC = temp.getDataC();
1188 } else {
1189 nodeDataC = arg.getDataC();
1190 }
1191 } else {
1192 nodeDataC = arg.getDataC();
1193 }
1194 switch(grad.getFunctionSpace().getTypeCode()) {
1195 case(Nodes):
1196 throw DudleyAdapterException("Error - Gradient at nodes is not supported.");
1197 break;
1198 case(ReducedNodes):
1199 throw DudleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1200 break;
1201 case(Elements):
1202 Dudley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1203 break;
1204 case(ReducedElements):
1205 Dudley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1206 break;
1207 case(FaceElements):
1208 Dudley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1209 break;
1210 case(ReducedFaceElements):
1211 Dudley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1212 break;
1213 case(Points):
1214 throw DudleyAdapterException("Error - Gradient at points is not supported.");
1215 break;
1216 case(DegreesOfFreedom):
1217 throw DudleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1218 break;
1219 case(ReducedDegreesOfFreedom):
1220 throw DudleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1221 break;
1222 default:
1223 stringstream temp;
1224 temp << "Error - Gradient: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1225 throw DudleyAdapterException(temp.str());
1226 break;
1227 }
1228 checkDudleyError();
1229 }
1230
1231 //
1232 // returns the size of elements:
1233 //
1234 void MeshAdapter::setToSize(escript::Data& size) const
1235 {
1236 Dudley_Mesh* mesh=m_dudleyMesh.get();
1237 escriptDataC tmp=size.getDataC();
1238 switch(size.getFunctionSpace().getTypeCode()) {
1239 case(Nodes):
1240 throw DudleyAdapterException("Error - Size of nodes is not supported.");
1241 break;
1242 case(ReducedNodes):
1243 throw DudleyAdapterException("Error - Size of reduced nodes is not supported.");
1244 break;
1245 case(Elements):
1246 Dudley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1247 break;
1248 case(ReducedElements):
1249 Dudley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1250 break;
1251 case(FaceElements):
1252 Dudley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1253 break;
1254 case(ReducedFaceElements):
1255 Dudley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1256 break;
1257 case(Points):
1258 throw DudleyAdapterException("Error - Size of point elements is not supported.");
1259 break;
1260 case(DegreesOfFreedom):
1261 throw DudleyAdapterException("Error - Size of degrees of freedom is not supported.");
1262 break;
1263 case(ReducedDegreesOfFreedom):
1264 throw DudleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1265 break;
1266 default:
1267 stringstream temp;
1268 temp << "Error - Element size: Dudley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1269 throw DudleyAdapterException(temp.str());
1270 break;
1271 }
1272 checkDudleyError();
1273 }
1274
1275 //
1276 // sets the location of nodes
1277 //
1278 void MeshAdapter::setNewX(const escript::Data& new_x)
1279 {
1280 Dudley_Mesh* mesh=m_dudleyMesh.get();
1281 escriptDataC tmp;
1282 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1283 if (newDomain!=*this)
1284 throw DudleyAdapterException("Error - Illegal domain of new point locations");
1285 if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {
1286 tmp = new_x.getDataC();
1287 Dudley_Mesh_setCoordinates(mesh,&tmp);
1288 } else {
1289 escript::Data new_x_inter=escript::Data( new_x, continuousFunction(asAbstractContinuousDomain()) );
1290 tmp = new_x_inter.getDataC();
1291 Dudley_Mesh_setCoordinates(mesh,&tmp);
1292 }
1293 checkDudleyError();
1294 }
1295
1296 //
1297 // Helper for the save* methods. Extracts optional data variable names and the
1298 // corresponding pointers from python dictionary. Caller must free arrays.
1299 //
1300 void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1301 {
1302 numData = boost::python::extract<int>(arg.attr("__len__")());
1303 /* win32 refactor */
1304 names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1305 data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1306 dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1307
1308 boost::python::list keys=arg.keys();
1309 for (int i=0; i<numData; ++i) {
1310 string n=boost::python::extract<string>(keys[i]);
1311 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1312 if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1313 throw DudleyAdapterException("Error: Data must be defined on same Domain");
1314 data[i] = d.getDataC();
1315 dataPtr[i] = &(data[i]);
1316 names[i] = TMPMEMALLOC(n.length()+1, char);
1317 strcpy(names[i], n.c_str());
1318 }
1319 }
1320
1321 //
1322 // saves mesh and optionally data arrays in openDX format
1323 //
1324 void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1325 {
1326 int num_data;
1327 char **names;
1328 escriptDataC *data;
1329 escriptDataC **ptr_data;
1330
1331 extractArgsFromDict(arg, num_data, names, data, ptr_data);
1332 Dudley_Mesh_saveDX(filename.c_str(), m_dudleyMesh.get(), num_data, names, ptr_data);
1333 checkDudleyError();
1334
1335 /* win32 refactor */
1336 TMPMEMFREE(data);
1337 TMPMEMFREE(ptr_data);
1338 for(int i=0; i<num_data; i++)
1339 {
1340 TMPMEMFREE(names[i]);
1341 }
1342 TMPMEMFREE(names);
1343
1344 return;
1345 }
1346
1347 //
1348 // saves mesh and optionally data arrays in VTK format
1349 //
1350 void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg, const string& metadata, const string& metadata_schema) const
1351 {
1352 int num_data;
1353 char **names;
1354 escriptDataC *data;
1355 escriptDataC **ptr_data;
1356
1357 extractArgsFromDict(arg, num_data, names, data, ptr_data);
1358 Dudley_Mesh_saveVTK(filename.c_str(), m_dudleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1359 checkDudleyError();
1360
1361 /* win32 refactor */
1362 TMPMEMFREE(data);
1363 TMPMEMFREE(ptr_data);
1364 for(int i=0; i<num_data; i++)
1365 {
1366 TMPMEMFREE(names[i]);
1367 }
1368 TMPMEMFREE(names);
1369 }
1370
1371 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1372 {
1373 #ifdef PASO_MPI
1374 index_t myFirstNode=0, myLastNode=0, k=0;
1375 index_t* globalNodeIndex=0;
1376 Dudley_Mesh* mesh_p=m_dudleyMesh.get();
1377 if (fs_code == DUDLEY_REDUCED_NODES)
1378 {
1379 myFirstNode = Dudley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1380 myLastNode = Dudley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1381 globalNodeIndex = Dudley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1382 }
1383 else
1384 {
1385 myFirstNode = Dudley_NodeFile_getFirstNode(mesh_p->Nodes);
1386 myLastNode = Dudley_NodeFile_getLastNode(mesh_p->Nodes);
1387 globalNodeIndex = Dudley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1388 }
1389 k=globalNodeIndex[id];
1390 return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1391 #endif
1392 return true;
1393 }
1394
1395
1396
1397 //
1398 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1399 //
1400 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1401 const int row_blocksize,
1402 const escript::FunctionSpace& row_functionspace,
1403 const int column_blocksize,
1404 const escript::FunctionSpace& column_functionspace,
1405 const int type) const
1406 {
1407 int reduceRowOrder=0;
1408 int reduceColOrder=0;
1409 // is the domain right?
1410 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1411 if (row_domain!=*this)
1412 throw DudleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1413 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1414 if (col_domain!=*this)
1415 throw DudleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1416 // is the function space type right
1417 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1418 reduceRowOrder=0;
1419 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1420 reduceRowOrder=1;
1421 } else {
1422 throw DudleyAdapterException("Error - illegal function space type for system matrix rows.");
1423 }
1424 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1425 reduceColOrder=0;
1426 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1427 reduceColOrder=1;
1428 } else {
1429 throw DudleyAdapterException("Error - illegal function space type for system matrix columns.");
1430 }
1431 // generate matrix:
1432
1433 Paso_SystemMatrixPattern* fsystemMatrixPattern=Dudley_getPattern(getDudley_Mesh(),reduceRowOrder,reduceColOrder);
1434 checkDudleyError();
1435 Paso_SystemMatrix* fsystemMatrix;
1436 int trilinos = 0;
1437 if (trilinos) {
1438 #ifdef TRILINOS
1439 /* Allocation Epetra_VrbMatrix here */
1440 #endif
1441 }
1442 else {
1443 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1444 }
1445 checkPasoError();
1446 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1447 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1448 }
1449
1450 //
1451 // creates a TransportProblemAdapter
1452 //
1453 TransportProblemAdapter MeshAdapter::newTransportProblem(
1454 const bool useBackwardEuler,
1455 const int blocksize,
1456 const escript::FunctionSpace& functionspace,
1457 const int type) const
1458 {
1459 int reduceOrder=0;
1460 // is the domain right?
1461 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1462 if (domain!=*this)
1463 throw DudleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1464 // is the function space type right
1465 if (functionspace.getTypeCode()==DegreesOfFreedom) {
1466 reduceOrder=0;
1467 } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1468 reduceOrder=1;
1469 } else {
1470 throw DudleyAdapterException("Error - illegal function space type for system matrix rows.");
1471 }
1472 // generate matrix:
1473
1474 Paso_SystemMatrixPattern* fsystemMatrixPattern=Dudley_getPattern(getDudley_Mesh(),reduceOrder,reduceOrder);
1475 checkDudleyError();
1476 Paso_TransportProblem* transportProblem;
1477 transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1478 checkPasoError();
1479 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1480 return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1481 }
1482
1483 //
1484 // vtkObject MeshAdapter::createVtkObject() const
1485 // TODO:
1486 //
1487
1488 //
1489 // returns true if data at the atom_type is considered as being cell centered:
1490 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1491 {
1492 switch(functionSpaceCode) {
1493 case(Nodes):
1494 case(DegreesOfFreedom):
1495 case(ReducedDegreesOfFreedom):
1496 return false;
1497 break;
1498 case(Elements):
1499 case(FaceElements):
1500 case(Points):
1501 case(ReducedElements):
1502 case(ReducedFaceElements):
1503 return true;
1504 break;
1505 default:
1506 stringstream temp;
1507 temp << "Error - Cell: Dudley does not know anything about function space type " << functionSpaceCode;
1508 throw DudleyAdapterException(temp.str());
1509 break;
1510 }
1511 checkDudleyError();
1512 return false;
1513 }
1514
1515 bool
1516 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1517 {
1518 /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1519 class 1: DOF <-> Nodes
1520 class 2: ReducedDOF <-> ReducedNodes
1521 class 3: Points
1522 class 4: Elements
1523 class 5: ReducedElements
1524 class 6: FaceElements
1525 class 7: ReducedFaceElements
1526 class 8: ContactElementZero <-> ContactElementOne
1527 class 9: ReducedContactElementZero <-> ReducedContactElementOne
1528
1529 There is also a set of lines. Interpolation is possible down a line but not between lines.
1530 class 1 and 2 belong to all lines so aren't considered.
1531 line 0: class 3
1532 line 1: class 4,5
1533 line 2: class 6,7
1534 line 3: class 8,9
1535
1536 For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1537 eg hasnodes is true if we have at least one instance of Nodes.
1538 */
1539 if (fs.empty())
1540 {
1541 return false;
1542 }
1543 vector<int> hasclass(10);
1544 vector<int> hasline(4);
1545 bool hasnodes=false;
1546 bool hasrednodes=false;
1547 for (int i=0;i<fs.size();++i)
1548 {
1549 switch(fs[i])
1550 {
1551 case(Nodes): hasnodes=true; // no break is deliberate
1552 case(DegreesOfFreedom):
1553 hasclass[1]=1;
1554 break;
1555 case(ReducedNodes): hasrednodes=true; // no break is deliberate
1556 case(ReducedDegreesOfFreedom):
1557 hasclass[2]=1;
1558 break;
1559 case(Points):
1560 hasline[0]=1;
1561 hasclass[3]=1;
1562 break;
1563 case(Elements):
1564 hasclass[4]=1;
1565 hasline[1]=1;
1566 break;
1567 case(ReducedElements):
1568 hasclass[5]=1;
1569 hasline[1]=1;
1570 break;
1571 case(FaceElements):
1572 hasclass[6]=1;
1573 hasline[2]=1;
1574 break;
1575 case(ReducedFaceElements):
1576 hasclass[7]=1;
1577 hasline[2]=1;
1578 break;
1579 default:
1580 return false;
1581 }
1582 }
1583 int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1584 // fail if we have more than one leaf group
1585
1586 if (totlines>1)
1587 {
1588 return false; // there are at least two branches we can't interpolate between
1589 }
1590 else if (totlines==1)
1591 {
1592 if (hasline[0]==1) // we have points
1593 {
1594 resultcode=Points;
1595 }
1596 else if (hasline[1]==1)
1597 {
1598 if (hasclass[5]==1)
1599 {
1600 resultcode=ReducedElements;
1601 }
1602 else
1603 {
1604 resultcode=Elements;
1605 }
1606 }
1607 else if (hasline[2]==1)
1608 {
1609 if (hasclass[7]==1)
1610 {
1611 resultcode=ReducedFaceElements;
1612 }
1613 else
1614 {
1615 resultcode=FaceElements;
1616 }
1617 }
1618 else // so we must be in line3
1619 {
1620
1621 throw DudleyAdapterException("Programmer Error - choosing between contact elements - we should never get here.");
1622
1623 }
1624 }
1625 else // totlines==0
1626 {
1627 if (hasclass[2]==1)
1628 {
1629 // something from class 2
1630 resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1631 }
1632 else
1633 { // something from class 1
1634 resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1635 }
1636 }
1637 return true;
1638 }
1639
1640 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1641 {
1642 switch(functionSpaceType_source) {
1643 case(Nodes):
1644 switch(functionSpaceType_target) {
1645 case(Nodes):
1646 case(ReducedNodes):
1647 case(ReducedDegreesOfFreedom):
1648 case(DegreesOfFreedom):
1649 case(Elements):
1650 case(ReducedElements):
1651 case(FaceElements):
1652 case(ReducedFaceElements):
1653 case(Points):
1654 return true;
1655 default:
1656 stringstream temp;
1657 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1658 throw DudleyAdapterException(temp.str());
1659 }
1660 break;
1661 case(ReducedNodes):
1662 switch(functionSpaceType_target) {
1663 case(ReducedNodes):
1664 case(ReducedDegreesOfFreedom):
1665 case(Elements):
1666 case(ReducedElements):
1667 case(FaceElements):
1668 case(ReducedFaceElements):
1669 case(Points):
1670 return true;
1671 case(Nodes):
1672 case(DegreesOfFreedom):
1673 return false;
1674 default:
1675 stringstream temp;
1676 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1677 throw DudleyAdapterException(temp.str());
1678 }
1679 break;
1680 case(Elements):
1681 if (functionSpaceType_target==Elements) {
1682 return true;
1683 } else if (functionSpaceType_target==ReducedElements) {
1684 return true;
1685 } else {
1686 return false;
1687 }
1688 case(ReducedElements):
1689 if (functionSpaceType_target==ReducedElements) {
1690 return true;
1691 } else {
1692 return false;
1693 }
1694 case(FaceElements):
1695 if (functionSpaceType_target==FaceElements) {
1696 return true;
1697 } else if (functionSpaceType_target==ReducedFaceElements) {
1698 return true;
1699 } else {
1700 return false;
1701 }
1702 case(ReducedFaceElements):
1703 if (functionSpaceType_target==ReducedFaceElements) {
1704 return true;
1705 } else {
1706 return false;
1707 }
1708 case(Points):
1709 if (functionSpaceType_target==Points) {
1710 return true;
1711 } else {
1712 return false;
1713 }
1714 case(DegreesOfFreedom):
1715 switch(functionSpaceType_target) {
1716 case(ReducedDegreesOfFreedom):
1717 case(DegreesOfFreedom):
1718 case(Nodes):
1719 case(ReducedNodes):
1720 case(Elements):
1721 case(ReducedElements):
1722 case(Points):
1723 case(FaceElements):
1724 case(ReducedFaceElements):
1725 return true;
1726 default:
1727 stringstream temp;
1728 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1729 throw DudleyAdapterException(temp.str());
1730 }
1731 break;
1732 case(ReducedDegreesOfFreedom):
1733 switch(functionSpaceType_target) {
1734 case(ReducedDegreesOfFreedom):
1735 case(ReducedNodes):
1736 case(Elements):
1737 case(ReducedElements):
1738 case(FaceElements):
1739 case(ReducedFaceElements):
1740 case(Points):
1741 return true;
1742 case(Nodes):
1743 case(DegreesOfFreedom):
1744 return false;
1745 default:
1746 stringstream temp;
1747 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1748 throw DudleyAdapterException(temp.str());
1749 }
1750 break;
1751 default:
1752 stringstream temp;
1753 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_source;
1754 throw DudleyAdapterException(temp.str());
1755 break;
1756 }
1757 checkDudleyError();
1758 return false;
1759 }
1760
1761 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1762 {
1763 return false;
1764 }
1765
1766 bool MeshAdapter::operator==(const AbstractDomain& other) const
1767 {
1768 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1769 if (temp!=0) {
1770 return (m_dudleyMesh==temp->m_dudleyMesh);
1771 } else {
1772 return false;
1773 }
1774 }
1775
1776 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1777 {
1778 return !(operator==(other));
1779 }
1780
1781 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1782 {
1783 Dudley_Mesh* mesh=m_dudleyMesh.get();
1784 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1785 checkPasoError();
1786 return out;
1787 }
1788 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1789 {
1790 Dudley_Mesh* mesh=m_dudleyMesh.get();
1791 int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1792 checkPasoError();
1793 return out;
1794 }
1795
1796 escript::Data MeshAdapter::getX() const
1797 {
1798 return continuousFunction(asAbstractContinuousDomain()).getX();
1799 }
1800
1801 escript::Data MeshAdapter::getNormal() const
1802 {
1803 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1804 }
1805
1806 escript::Data MeshAdapter::getSize() const
1807 {
1808 return escript::function(asAbstractContinuousDomain()).getSize();
1809 }
1810
1811 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1812 {
1813 int *out = NULL;
1814 Dudley_Mesh* mesh=m_dudleyMesh.get();
1815 switch (functionSpaceType) {
1816 case(Nodes):
1817 out=mesh->Nodes->Id;
1818 break;
1819 case(ReducedNodes):
1820 out=mesh->Nodes->reducedNodesId;
1821 break;
1822 case(Elements):
1823 out=mesh->Elements->Id;
1824 break;
1825 case(ReducedElements):
1826 out=mesh->Elements->Id;
1827 break;
1828 case(FaceElements):
1829 out=mesh->FaceElements->Id;
1830 break;
1831 case(ReducedFaceElements):
1832 out=mesh->FaceElements->Id;
1833 break;
1834 case(Points):
1835 out=mesh->Points->Id;
1836 break;
1837 case(DegreesOfFreedom):
1838 out=mesh->Nodes->degreesOfFreedomId;
1839 break;
1840 case(ReducedDegreesOfFreedom):
1841 out=mesh->Nodes->reducedDegreesOfFreedomId;
1842 break;
1843 default:
1844 stringstream temp;
1845 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1846 throw DudleyAdapterException(temp.str());
1847 break;
1848 }
1849 return out;
1850 }
1851 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1852 {
1853 int out=0;
1854 Dudley_Mesh* mesh=m_dudleyMesh.get();
1855 switch (functionSpaceType) {
1856 case(Nodes):
1857 out=mesh->Nodes->Tag[sampleNo];
1858 break;
1859 case(ReducedNodes):
1860 throw DudleyAdapterException(" Error - ReducedNodes does not support tags.");
1861 break;
1862 case(Elements):
1863 out=mesh->Elements->Tag[sampleNo];
1864 break;
1865 case(ReducedElements):
1866 out=mesh->Elements->Tag[sampleNo];
1867 break;
1868 case(FaceElements):
1869 out=mesh->FaceElements->Tag[sampleNo];
1870 break;
1871 case(ReducedFaceElements):
1872 out=mesh->FaceElements->Tag[sampleNo];
1873 break;
1874 case(Points):
1875 out=mesh->Points->Tag[sampleNo];
1876 break;
1877 case(DegreesOfFreedom):
1878 throw DudleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1879 break;
1880 case(ReducedDegreesOfFreedom):
1881 throw DudleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1882 break;
1883 default:
1884 stringstream temp;
1885 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1886 throw DudleyAdapterException(temp.str());
1887 break;
1888 }
1889 return out;
1890 }
1891
1892
1893 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1894 {
1895 Dudley_Mesh* mesh=m_dudleyMesh.get();
1896 escriptDataC tmp=mask.getDataC();
1897 switch(functionSpaceType) {
1898 case(Nodes):
1899 Dudley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1900 break;
1901 case(ReducedNodes):
1902 throw DudleyAdapterException("Error - ReducedNodes does not support tags");
1903 break;
1904 case(DegreesOfFreedom):
1905 throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
1906 break;
1907 case(ReducedDegreesOfFreedom):
1908 throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1909 break;
1910 case(Elements):
1911 Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1912 break;
1913 case(ReducedElements):
1914 Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1915 break;
1916 case(FaceElements):
1917 Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1918 break;
1919 case(ReducedFaceElements):
1920 Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1921 break;
1922 case(Points):
1923 Dudley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1924 break;
1925 default:
1926 stringstream temp;
1927 temp << "Error - Dudley does not know anything about function space type " << functionSpaceType;
1928 throw DudleyAdapterException(temp.str());
1929 }
1930 checkDudleyError();
1931 return;
1932 }
1933
1934 void MeshAdapter::setTagMap(const string& name, int tag)
1935 {
1936 Dudley_Mesh* mesh=m_dudleyMesh.get();
1937 Dudley_Mesh_addTagMap(mesh, name.c_str(),tag);
1938 checkPasoError();
1939 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1940 }
1941
1942 int MeshAdapter::getTag(const string& name) const
1943 {
1944 Dudley_Mesh* mesh=m_dudleyMesh.get();
1945 int tag=0;
1946 tag=Dudley_Mesh_getTag(mesh, name.c_str());
1947 checkPasoError();
1948 // throwStandardException("MeshAdapter::getTag is not implemented.");
1949 return tag;
1950 }
1951
1952 bool MeshAdapter::isValidTagName(const string& name) const
1953 {
1954 Dudley_Mesh* mesh=m_dudleyMesh.get();
1955 return Dudley_Mesh_isValidTagName(mesh,name.c_str());
1956 }
1957
1958 string MeshAdapter::showTagNames() const
1959 {
1960 stringstream temp;
1961 Dudley_Mesh* mesh=m_dudleyMesh.get();
1962 Dudley_TagMap* tag_map=mesh->TagMap;
1963 while (tag_map) {
1964 temp << tag_map->name;
1965 tag_map=tag_map->next;
1966 if (tag_map) temp << ", ";
1967 }
1968 return temp.str();
1969 }
1970
1971 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
1972 {
1973 Dudley_Mesh* mesh=m_dudleyMesh.get();
1974 dim_t numTags=0;
1975 switch(functionSpaceCode) {
1976 case(Nodes):
1977 numTags=mesh->Nodes->numTagsInUse;
1978 break;
1979 case(ReducedNodes):
1980 throw DudleyAdapterException("Error - ReducedNodes does not support tags");
1981 break;
1982 case(DegreesOfFreedom):
1983 throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
1984 break;
1985 case(ReducedDegreesOfFreedom):
1986 throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1987 break;
1988 case(Elements):
1989 case(ReducedElements):
1990 numTags=mesh->Elements->numTagsInUse;
1991 break;
1992 case(FaceElements):
1993 case(ReducedFaceElements):
1994 numTags=mesh->FaceElements->numTagsInUse;
1995 break;
1996 case(Points):
1997 numTags=mesh->Points->numTagsInUse;
1998 break;
1999 default:
2000 stringstream temp;
2001 temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2002 throw DudleyAdapterException(temp.str());
2003 }
2004 return numTags;
2005 }
2006
2007 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2008 {
2009 Dudley_Mesh* mesh=m_dudleyMesh.get();
2010 index_t* tags=NULL;
2011 switch(functionSpaceCode) {
2012 case(Nodes):
2013 tags=mesh->Nodes->tagsInUse;
2014 break;
2015 case(ReducedNodes):
2016 throw DudleyAdapterException("Error - ReducedNodes does not support tags");
2017 break;
2018 case(DegreesOfFreedom):
2019 throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
2020 break;
2021 case(ReducedDegreesOfFreedom):
2022 throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2023 break;
2024 case(Elements):
2025 case(ReducedElements):
2026 tags=mesh->Elements->tagsInUse;
2027 break;
2028 case(FaceElements):
2029 case(ReducedFaceElements):
2030 tags=mesh->FaceElements->tagsInUse;
2031 break;
2032 case(Points):
2033 tags=mesh->Points->tagsInUse;
2034 break;
2035 default:
2036 stringstream temp;
2037 temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2038 throw DudleyAdapterException(temp.str());
2039 }
2040 return tags;
2041 }
2042
2043
2044 bool MeshAdapter::canTag(int functionSpaceCode) const
2045 {
2046 switch(functionSpaceCode) {
2047 case(Nodes):
2048 case(Elements):
2049 case(ReducedElements):
2050 case(FaceElements):
2051 case(ReducedFaceElements):
2052 case(Points):
2053 return true;
2054 case(ReducedNodes):
2055 case(DegreesOfFreedom):
2056 case(ReducedDegreesOfFreedom):
2057 return false;
2058 default:
2059 return false;
2060 }
2061 }
2062
2063 AbstractDomain::StatusType MeshAdapter::getStatus() const
2064 {
2065 Dudley_Mesh* mesh=m_dudleyMesh.get();
2066 return Dudley_Mesh_getStatus(mesh);
2067 }
2068
2069 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2070 {
2071
2072 Dudley_Mesh* mesh=m_dudleyMesh.get();
2073 int order =-1;
2074 switch(functionSpaceCode) {
2075 case(Nodes):
2076 case(DegreesOfFreedom):
2077 order=mesh->approximationOrder;
2078 break;
2079 case(ReducedNodes):
2080 case(ReducedDegreesOfFreedom):
2081 order=mesh->reducedApproximationOrder;
2082 break;
2083 case(Elements):
2084 case(FaceElements):
2085 case(Points):
2086 order=mesh->integrationOrder;
2087 break;
2088 case(ReducedElements):
2089 case(ReducedFaceElements):
2090 order=mesh->reducedIntegrationOrder;
2091 break;
2092 default:
2093 stringstream temp;
2094 temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2095 throw DudleyAdapterException(temp.str());
2096 }
2097 return order;
2098 }
2099
2100
2101 bool MeshAdapter::supportsContactElements() const
2102 {
2103 return false;
2104 }
2105
2106 ReferenceElementSetWrapper::ReferenceElementSetWrapper(ElementTypeId id, index_t order, index_t reducedOrder)
2107 {
2108 m_refSet = Dudley_ReferenceElementSet_alloc(id, order, reducedOrder);
2109 }
2110
2111 ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2112 {
2113 Dudley_ReferenceElementSet_dealloc(m_refSet);
2114 }
2115
2116 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26