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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2989 - (hide annotations)
Thu Mar 18 01:45:55 2010 UTC (9 years ago) by gross
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 85811 byte(s)
- bug in error handling in the case of lumped mass matrix fixed.
- switched back to rowsum lumping as there are problems with HRZ lumping for order 2 


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