/[escript]/branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Contents of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4478 - (show annotations)
Thu Jun 20 06:34:02 2013 UTC (5 years, 8 months ago) by caltinay
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 91385 byte(s)
Moved appendRankToFilename to header within namespace to remove code duplication
and fixed memory leaks this way.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26