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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26