/[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 3239 - (show annotations)
Tue Oct 5 05:31:20 2010 UTC (8 years, 10 months ago) by jfenwick
File size: 72197 byte(s)
Put the tsunami example back to use finley.
Moved a bunch of method signatures up into the AbstractDomain class

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 ASM_ptr 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 SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace, column_blocksize,column_functionspace);
1445 return ASM_ptr(sma);
1446 }
1447
1448 //
1449 // creates a TransportProblemAdapter
1450 //
1451 ATP_ptr MeshAdapter::newTransportProblem(
1452 const bool useBackwardEuler,
1453 const int blocksize,
1454 const escript::FunctionSpace& functionspace,
1455 const int type) const
1456 {
1457 int reduceOrder=0;
1458 // is the domain right?
1459 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1460 if (domain!=*this)
1461 throw DudleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1462 // is the function space type right
1463 if (functionspace.getTypeCode()==DegreesOfFreedom) {
1464 reduceOrder=0;
1465 } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1466 reduceOrder=1;
1467 } else {
1468 throw DudleyAdapterException("Error - illegal function space type for system matrix rows.");
1469 }
1470 // generate matrix:
1471
1472 Paso_SystemMatrixPattern* fsystemMatrixPattern=Dudley_getPattern(getDudley_Mesh(),reduceOrder,reduceOrder);
1473 checkDudleyError();
1474 Paso_TransportProblem* transportProblem;
1475 transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1476 checkPasoError();
1477 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1478 AbstractTransportProblem* atp=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1479 return ATP_ptr(atp);
1480 }
1481
1482 //
1483 // vtkObject MeshAdapter::createVtkObject() const
1484 // TODO:
1485 //
1486
1487 //
1488 // returns true if data at the atom_type is considered as being cell centered:
1489 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1490 {
1491 switch(functionSpaceCode) {
1492 case(Nodes):
1493 case(DegreesOfFreedom):
1494 case(ReducedDegreesOfFreedom):
1495 return false;
1496 break;
1497 case(Elements):
1498 case(FaceElements):
1499 case(Points):
1500 case(ReducedElements):
1501 case(ReducedFaceElements):
1502 return true;
1503 break;
1504 default:
1505 stringstream temp;
1506 temp << "Error - Cell: Dudley does not know anything about function space type " << functionSpaceCode;
1507 throw DudleyAdapterException(temp.str());
1508 break;
1509 }
1510 checkDudleyError();
1511 return false;
1512 }
1513
1514 bool
1515 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1516 {
1517 /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1518 class 1: DOF <-> Nodes
1519 class 2: ReducedDOF <-> ReducedNodes
1520 class 3: Points
1521 class 4: Elements
1522 class 5: ReducedElements
1523 class 6: FaceElements
1524 class 7: ReducedFaceElements
1525 class 8: ContactElementZero <-> ContactElementOne
1526 class 9: ReducedContactElementZero <-> ReducedContactElementOne
1527
1528 There is also a set of lines. Interpolation is possible down a line but not between lines.
1529 class 1 and 2 belong to all lines so aren't considered.
1530 line 0: class 3
1531 line 1: class 4,5
1532 line 2: class 6,7
1533 line 3: class 8,9
1534
1535 For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1536 eg hasnodes is true if we have at least one instance of Nodes.
1537 */
1538 if (fs.empty())
1539 {
1540 return false;
1541 }
1542 vector<int> hasclass(10);
1543 vector<int> hasline(4);
1544 bool hasnodes=false;
1545 bool hasrednodes=false;
1546 for (int i=0;i<fs.size();++i)
1547 {
1548 switch(fs[i])
1549 {
1550 case(Nodes): hasnodes=true; // no break is deliberate
1551 case(DegreesOfFreedom):
1552 hasclass[1]=1;
1553 break;
1554 case(ReducedNodes): hasrednodes=true; // no break is deliberate
1555 case(ReducedDegreesOfFreedom):
1556 hasclass[2]=1;
1557 break;
1558 case(Points):
1559 hasline[0]=1;
1560 hasclass[3]=1;
1561 break;
1562 case(Elements):
1563 hasclass[4]=1;
1564 hasline[1]=1;
1565 break;
1566 case(ReducedElements):
1567 hasclass[5]=1;
1568 hasline[1]=1;
1569 break;
1570 case(FaceElements):
1571 hasclass[6]=1;
1572 hasline[2]=1;
1573 break;
1574 case(ReducedFaceElements):
1575 hasclass[7]=1;
1576 hasline[2]=1;
1577 break;
1578 default:
1579 return false;
1580 }
1581 }
1582 int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1583 // fail if we have more than one leaf group
1584
1585 if (totlines>1)
1586 {
1587 return false; // there are at least two branches we can't interpolate between
1588 }
1589 else if (totlines==1)
1590 {
1591 if (hasline[0]==1) // we have points
1592 {
1593 resultcode=Points;
1594 }
1595 else if (hasline[1]==1)
1596 {
1597 if (hasclass[5]==1)
1598 {
1599 resultcode=ReducedElements;
1600 }
1601 else
1602 {
1603 resultcode=Elements;
1604 }
1605 }
1606 else if (hasline[2]==1)
1607 {
1608 if (hasclass[7]==1)
1609 {
1610 resultcode=ReducedFaceElements;
1611 }
1612 else
1613 {
1614 resultcode=FaceElements;
1615 }
1616 }
1617 else // so we must be in line3
1618 {
1619
1620 throw DudleyAdapterException("Programmer Error - choosing between contact elements - we should never get here.");
1621
1622 }
1623 }
1624 else // totlines==0
1625 {
1626 if (hasclass[2]==1)
1627 {
1628 // something from class 2
1629 resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1630 }
1631 else
1632 { // something from class 1
1633 resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1634 }
1635 }
1636 return true;
1637 }
1638
1639 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1640 {
1641 switch(functionSpaceType_source) {
1642 case(Nodes):
1643 switch(functionSpaceType_target) {
1644 case(Nodes):
1645 case(ReducedNodes):
1646 case(ReducedDegreesOfFreedom):
1647 case(DegreesOfFreedom):
1648 case(Elements):
1649 case(ReducedElements):
1650 case(FaceElements):
1651 case(ReducedFaceElements):
1652 case(Points):
1653 return true;
1654 default:
1655 stringstream temp;
1656 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1657 throw DudleyAdapterException(temp.str());
1658 }
1659 break;
1660 case(ReducedNodes):
1661 switch(functionSpaceType_target) {
1662 case(ReducedNodes):
1663 case(ReducedDegreesOfFreedom):
1664 case(Elements):
1665 case(ReducedElements):
1666 case(FaceElements):
1667 case(ReducedFaceElements):
1668 case(Points):
1669 return true;
1670 case(Nodes):
1671 case(DegreesOfFreedom):
1672 return false;
1673 default:
1674 stringstream temp;
1675 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1676 throw DudleyAdapterException(temp.str());
1677 }
1678 break;
1679 case(Elements):
1680 if (functionSpaceType_target==Elements) {
1681 return true;
1682 } else if (functionSpaceType_target==ReducedElements) {
1683 return true;
1684 } else {
1685 return false;
1686 }
1687 case(ReducedElements):
1688 if (functionSpaceType_target==ReducedElements) {
1689 return true;
1690 } else {
1691 return false;
1692 }
1693 case(FaceElements):
1694 if (functionSpaceType_target==FaceElements) {
1695 return true;
1696 } else if (functionSpaceType_target==ReducedFaceElements) {
1697 return true;
1698 } else {
1699 return false;
1700 }
1701 case(ReducedFaceElements):
1702 if (functionSpaceType_target==ReducedFaceElements) {
1703 return true;
1704 } else {
1705 return false;
1706 }
1707 case(Points):
1708 if (functionSpaceType_target==Points) {
1709 return true;
1710 } else {
1711 return false;
1712 }
1713 case(DegreesOfFreedom):
1714 switch(functionSpaceType_target) {
1715 case(ReducedDegreesOfFreedom):
1716 case(DegreesOfFreedom):
1717 case(Nodes):
1718 case(ReducedNodes):
1719 case(Elements):
1720 case(ReducedElements):
1721 case(Points):
1722 case(FaceElements):
1723 case(ReducedFaceElements):
1724 return true;
1725 default:
1726 stringstream temp;
1727 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1728 throw DudleyAdapterException(temp.str());
1729 }
1730 break;
1731 case(ReducedDegreesOfFreedom):
1732 switch(functionSpaceType_target) {
1733 case(ReducedDegreesOfFreedom):
1734 case(ReducedNodes):
1735 case(Elements):
1736 case(ReducedElements):
1737 case(FaceElements):
1738 case(ReducedFaceElements):
1739 case(Points):
1740 return true;
1741 case(Nodes):
1742 case(DegreesOfFreedom):
1743 return false;
1744 default:
1745 stringstream temp;
1746 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1747 throw DudleyAdapterException(temp.str());
1748 }
1749 break;
1750 default:
1751 stringstream temp;
1752 temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_source;
1753 throw DudleyAdapterException(temp.str());
1754 break;
1755 }
1756 checkDudleyError();
1757 return false;
1758 }
1759
1760 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1761 {
1762 return false;
1763 }
1764
1765 bool MeshAdapter::operator==(const AbstractDomain& other) const
1766 {
1767 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1768 if (temp!=0) {
1769 return (m_dudleyMesh==temp->m_dudleyMesh);
1770 } else {
1771 return false;
1772 }
1773 }
1774
1775 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1776 {
1777 return !(operator==(other));
1778 }
1779
1780 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1781 {
1782 Dudley_Mesh* mesh=m_dudleyMesh.get();
1783 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1784 checkPasoError();
1785 return out;
1786 }
1787 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1788 {
1789 Dudley_Mesh* mesh=m_dudleyMesh.get();
1790 int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1791 checkPasoError();
1792 return out;
1793 }
1794
1795 escript::Data MeshAdapter::getX() const
1796 {
1797 return continuousFunction(asAbstractContinuousDomain()).getX();
1798 }
1799
1800 escript::Data MeshAdapter::getNormal() const
1801 {
1802 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1803 }
1804
1805 escript::Data MeshAdapter::getSize() const
1806 {
1807 return escript::function(asAbstractContinuousDomain()).getSize();
1808 }
1809
1810 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1811 {
1812 int *out = NULL;
1813 Dudley_Mesh* mesh=m_dudleyMesh.get();
1814 switch (functionSpaceType) {
1815 case(Nodes):
1816 out=mesh->Nodes->Id;
1817 break;
1818 case(ReducedNodes):
1819 out=mesh->Nodes->reducedNodesId;
1820 break;
1821 case(Elements):
1822 out=mesh->Elements->Id;
1823 break;
1824 case(ReducedElements):
1825 out=mesh->Elements->Id;
1826 break;
1827 case(FaceElements):
1828 out=mesh->FaceElements->Id;
1829 break;
1830 case(ReducedFaceElements):
1831 out=mesh->FaceElements->Id;
1832 break;
1833 case(Points):
1834 out=mesh->Points->Id;
1835 break;
1836 case(DegreesOfFreedom):
1837 out=mesh->Nodes->degreesOfFreedomId;
1838 break;
1839 case(ReducedDegreesOfFreedom):
1840 out=mesh->Nodes->reducedDegreesOfFreedomId;
1841 break;
1842 default:
1843 stringstream temp;
1844 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1845 throw DudleyAdapterException(temp.str());
1846 break;
1847 }
1848 return out;
1849 }
1850 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1851 {
1852 int out=0;
1853 Dudley_Mesh* mesh=m_dudleyMesh.get();
1854 switch (functionSpaceType) {
1855 case(Nodes):
1856 out=mesh->Nodes->Tag[sampleNo];
1857 break;
1858 case(ReducedNodes):
1859 throw DudleyAdapterException(" Error - ReducedNodes does not support tags.");
1860 break;
1861 case(Elements):
1862 out=mesh->Elements->Tag[sampleNo];
1863 break;
1864 case(ReducedElements):
1865 out=mesh->Elements->Tag[sampleNo];
1866 break;
1867 case(FaceElements):
1868 out=mesh->FaceElements->Tag[sampleNo];
1869 break;
1870 case(ReducedFaceElements):
1871 out=mesh->FaceElements->Tag[sampleNo];
1872 break;
1873 case(Points):
1874 out=mesh->Points->Tag[sampleNo];
1875 break;
1876 case(DegreesOfFreedom):
1877 throw DudleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1878 break;
1879 case(ReducedDegreesOfFreedom):
1880 throw DudleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1881 break;
1882 default:
1883 stringstream temp;
1884 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1885 throw DudleyAdapterException(temp.str());
1886 break;
1887 }
1888 return out;
1889 }
1890
1891
1892 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1893 {
1894 Dudley_Mesh* mesh=m_dudleyMesh.get();
1895 escriptDataC tmp=mask.getDataC();
1896 switch(functionSpaceType) {
1897 case(Nodes):
1898 Dudley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1899 break;
1900 case(ReducedNodes):
1901 throw DudleyAdapterException("Error - ReducedNodes does not support tags");
1902 break;
1903 case(DegreesOfFreedom):
1904 throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
1905 break;
1906 case(ReducedDegreesOfFreedom):
1907 throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1908 break;
1909 case(Elements):
1910 Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1911 break;
1912 case(ReducedElements):
1913 Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1914 break;
1915 case(FaceElements):
1916 Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1917 break;
1918 case(ReducedFaceElements):
1919 Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1920 break;
1921 case(Points):
1922 Dudley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1923 break;
1924 default:
1925 stringstream temp;
1926 temp << "Error - Dudley does not know anything about function space type " << functionSpaceType;
1927 throw DudleyAdapterException(temp.str());
1928 }
1929 checkDudleyError();
1930 return;
1931 }
1932
1933 void MeshAdapter::setTagMap(const string& name, int tag)
1934 {
1935 Dudley_Mesh* mesh=m_dudleyMesh.get();
1936 Dudley_Mesh_addTagMap(mesh, name.c_str(),tag);
1937 checkPasoError();
1938 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1939 }
1940
1941 int MeshAdapter::getTag(const string& name) const
1942 {
1943 Dudley_Mesh* mesh=m_dudleyMesh.get();
1944 int tag=0;
1945 tag=Dudley_Mesh_getTag(mesh, name.c_str());
1946 checkPasoError();
1947 // throwStandardException("MeshAdapter::getTag is not implemented.");
1948 return tag;
1949 }
1950
1951 bool MeshAdapter::isValidTagName(const string& name) const
1952 {
1953 Dudley_Mesh* mesh=m_dudleyMesh.get();
1954 return Dudley_Mesh_isValidTagName(mesh,name.c_str());
1955 }
1956
1957 string MeshAdapter::showTagNames() const
1958 {
1959 stringstream temp;
1960 Dudley_Mesh* mesh=m_dudleyMesh.get();
1961 Dudley_TagMap* tag_map=mesh->TagMap;
1962 while (tag_map) {
1963 temp << tag_map->name;
1964 tag_map=tag_map->next;
1965 if (tag_map) temp << ", ";
1966 }
1967 return temp.str();
1968 }
1969
1970 int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
1971 {
1972 Dudley_Mesh* mesh=m_dudleyMesh.get();
1973 dim_t numTags=0;
1974 switch(functionSpaceCode) {
1975 case(Nodes):
1976 numTags=mesh->Nodes->numTagsInUse;
1977 break;
1978 case(ReducedNodes):
1979 throw DudleyAdapterException("Error - ReducedNodes does not support tags");
1980 break;
1981 case(DegreesOfFreedom):
1982 throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
1983 break;
1984 case(ReducedDegreesOfFreedom):
1985 throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1986 break;
1987 case(Elements):
1988 case(ReducedElements):
1989 numTags=mesh->Elements->numTagsInUse;
1990 break;
1991 case(FaceElements):
1992 case(ReducedFaceElements):
1993 numTags=mesh->FaceElements->numTagsInUse;
1994 break;
1995 case(Points):
1996 numTags=mesh->Points->numTagsInUse;
1997 break;
1998 default:
1999 stringstream temp;
2000 temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2001 throw DudleyAdapterException(temp.str());
2002 }
2003 return numTags;
2004 }
2005
2006 const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2007 {
2008 Dudley_Mesh* mesh=m_dudleyMesh.get();
2009 index_t* tags=NULL;
2010 switch(functionSpaceCode) {
2011 case(Nodes):
2012 tags=mesh->Nodes->tagsInUse;
2013 break;
2014 case(ReducedNodes):
2015 throw DudleyAdapterException("Error - ReducedNodes does not support tags");
2016 break;
2017 case(DegreesOfFreedom):
2018 throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
2019 break;
2020 case(ReducedDegreesOfFreedom):
2021 throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2022 break;
2023 case(Elements):
2024 case(ReducedElements):
2025 tags=mesh->Elements->tagsInUse;
2026 break;
2027 case(FaceElements):
2028 case(ReducedFaceElements):
2029 tags=mesh->FaceElements->tagsInUse;
2030 break;
2031 case(Points):
2032 tags=mesh->Points->tagsInUse;
2033 break;
2034 default:
2035 stringstream temp;
2036 temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2037 throw DudleyAdapterException(temp.str());
2038 }
2039 return tags;
2040 }
2041
2042
2043 bool MeshAdapter::canTag(int functionSpaceCode) const
2044 {
2045 switch(functionSpaceCode) {
2046 case(Nodes):
2047 case(Elements):
2048 case(ReducedElements):
2049 case(FaceElements):
2050 case(ReducedFaceElements):
2051 case(Points):
2052 return true;
2053 case(ReducedNodes):
2054 case(DegreesOfFreedom):
2055 case(ReducedDegreesOfFreedom):
2056 return false;
2057 default:
2058 return false;
2059 }
2060 }
2061
2062 AbstractDomain::StatusType MeshAdapter::getStatus() const
2063 {
2064 Dudley_Mesh* mesh=m_dudleyMesh.get();
2065 return Dudley_Mesh_getStatus(mesh);
2066 }
2067
2068 int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2069 {
2070
2071 Dudley_Mesh* mesh=m_dudleyMesh.get();
2072 int order =-1;
2073 switch(functionSpaceCode) {
2074 case(Nodes):
2075 case(DegreesOfFreedom):
2076 order=mesh->approximationOrder;
2077 break;
2078 case(ReducedNodes):
2079 case(ReducedDegreesOfFreedom):
2080 order=mesh->reducedApproximationOrder;
2081 break;
2082 case(Elements):
2083 case(FaceElements):
2084 case(Points):
2085 order=mesh->integrationOrder;
2086 break;
2087 case(ReducedElements):
2088 case(ReducedFaceElements):
2089 order=mesh->reducedIntegrationOrder;
2090 break;
2091 default:
2092 stringstream temp;
2093 temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
2094 throw DudleyAdapterException(temp.str());
2095 }
2096 return order;
2097 }
2098
2099
2100 bool MeshAdapter::supportsContactElements() const
2101 {
2102 return false;
2103 }
2104
2105 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26