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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26