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

Contents of /branches/symbolic_from_3470/dudley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3789 - (show annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 2 months ago) by caltinay
File size: 72772 byte(s)
Fast forward to latest trunk revision 3788.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26