/[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 3269 - (show annotations)
Wed Oct 13 03:21:50 2010 UTC (10 years, 6 months ago) by jfenwick
File size: 72761 byte(s)
Fixed some intel compiler warnings.
Put linearPDEs back the way it was and present a common interface for dudley and finley (as per Lutz)

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26