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

Annotation of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3998 - (hide annotations)
Thu Sep 27 01:17:28 2012 UTC (6 years, 8 months ago) by caltinay
File size: 91311 byte(s)
Finley now throws in ownSample() for reduced function with mpisize>1.

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