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

Contents of /trunk/dudley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4114 - (show annotations)
Fri Dec 14 04:24:46 2012 UTC (6 years, 8 months ago) by caltinay
File size: 71448 byte(s)
Time to remove deprecated saveVTK/DX methods from Data and Domain.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26