/[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 2898 - (hide annotations)
Mon Feb 1 01:23:46 2010 UTC (9 years, 1 month ago) by caltinay
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 85765 byte(s)
Reverted the move of system_dep.h since it didn't solve the problem nicely.
Created a wrapper class for the needed functionality instead.

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     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
815    
816     checkFinleyError();
817     }
818    
819    
820 jgs 82 //
821 jgs 102 // adds linear PDE of second order into the right hand side only
822     //
823 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
824 jgs 102 {
825     Finley_Mesh* mesh=m_finleyMesh.get();
826 jgs 147
827 gross 798 escriptDataC _rhs=rhs.getDataC();
828     escriptDataC _X=X.getDataC();
829     escriptDataC _Y=Y.getDataC();
830     escriptDataC _y=y.getDataC();
831     escriptDataC _y_contact=y_contact.getDataC();
832    
833     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
834 jgs 102 checkFinleyError();
835 jgs 147
836 gross 798 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
837     checkFinleyError();
838 jgs 147
839 gross 798 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
840 jgs 102 checkFinleyError();
841     }
842 gross 1367 //
843     // adds PDE of second order into a transport problem
844     //
845     void MeshAdapter::addPDEToTransportProblem(
846 phornby 1628 TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,
847     const escript::Data& A, const escript::Data& B, const escript::Data& C,
848     const escript::Data& D,const escript::Data& X,const escript::Data& Y,
849     const escript::Data& d, const escript::Data& y,
850     const escript::Data& d_contact,const escript::Data& y_contact) const
851 gross 1367 {
852 jfenwick 1796 DataTypes::ShapeType shape;
853 gross 1370 source.expand();
854 gross 1367 escriptDataC _source=source.getDataC();
855     escriptDataC _M=M.getDataC();
856     escriptDataC _A=A.getDataC();
857     escriptDataC _B=B.getDataC();
858     escriptDataC _C=C.getDataC();
859     escriptDataC _D=D.getDataC();
860     escriptDataC _X=X.getDataC();
861     escriptDataC _Y=Y.getDataC();
862     escriptDataC _d=d.getDataC();
863     escriptDataC _y=y.getDataC();
864     escriptDataC _d_contact=d_contact.getDataC();
865     escriptDataC _y_contact=y_contact.getDataC();
866 jgs 149
867 gross 1367 Finley_Mesh* mesh=m_finleyMesh.get();
868     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();
869 phornby 1628
870 gross 1407 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
871 gross 1367 checkFinleyError();
872    
873 gross 1407 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
874 gross 1367 checkFinleyError();
875    
876     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
877     checkFinleyError();
878    
879     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
880     checkFinleyError();
881     }
882    
883 jgs 102 //
884 jgs 82 // interpolates data between different function spaces:
885     //
886 jgs 480 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
887 jgs 82 {
888 jfenwick 1872 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
889     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
890 phornby 1628 if (inDomain!=*this)
891     throw FinleyAdapterException("Error - Illegal domain of interpolant.");
892     if (targetDomain!=*this)
893     throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
894 jgs 82
895 phornby 1628 Finley_Mesh* mesh=m_finleyMesh.get();
896     escriptDataC _target=target.getDataC();
897     escriptDataC _in=in.getDataC();
898     switch(in.getFunctionSpace().getTypeCode()) {
899     case(Nodes):
900 gross 2385 switch(target.getFunctionSpace().getTypeCode()) {
901     case(Nodes):
902     case(ReducedNodes):
903     case(DegreesOfFreedom):
904     case(ReducedDegreesOfFreedom):
905     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
906 phornby 1628 break;
907 gross 2385 case(Elements):
908     case(ReducedElements):
909     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
910 phornby 1628 break;
911 gross 2385 case(FaceElements):
912     case(ReducedFaceElements):
913     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
914     break;
915     case(Points):
916     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
917     break;
918     case(ContactElementsZero):
919     case(ReducedContactElementsZero):
920     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
921     break;
922     case(ContactElementsOne):
923     case(ReducedContactElementsOne):
924     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
925     break;
926     default:
927     stringstream temp;
928     temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
929     throw FinleyAdapterException(temp.str());
930     break;
931     }
932     break;
933 phornby 1628 case(ReducedNodes):
934 gross 2385 switch(target.getFunctionSpace().getTypeCode()) {
935     case(Nodes):
936     case(ReducedNodes):
937     case(DegreesOfFreedom):
938     case(ReducedDegreesOfFreedom):
939 phornby 1628 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
940 gross 2385 break;
941     case(Elements):
942     case(ReducedElements):
943 phornby 1628 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
944 gross 2385 break;
945     case(FaceElements):
946     case(ReducedFaceElements):
947 phornby 1628 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
948 gross 2385 break;
949     case(Points):
950 phornby 1628 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
951 gross 2385 break;
952     case(ContactElementsZero):
953     case(ReducedContactElementsZero):
954 phornby 1628 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
955     break;
956 gross 2385 case(ContactElementsOne):
957     case(ReducedContactElementsOne):
958     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
959     break;
960     default:
961     stringstream temp;
962     temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
963     throw FinleyAdapterException(temp.str());
964     break;
965     }
966     break;
967 phornby 1628 case(Elements):
968 gross 2385 if (target.getFunctionSpace().getTypeCode()==Elements) {
969     Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
970     } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
971     Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
972     } else {
973     throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
974     }
975     break;
976 phornby 1628 case(ReducedElements):
977 gross 2385 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
978     Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
979     } else {
980     throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
981     }
982     break;
983 phornby 1628 case(FaceElements):
984 gross 2385 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
985     Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
986     } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
987     Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
988     } else {
989     throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
990     }
991     break;
992 phornby 1628 case(ReducedFaceElements):
993 gross 2385 if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
994     Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
995     } else {
996     throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
997     }
998     break;
999 phornby 1628 case(Points):
1000 gross 2385 if (target.getFunctionSpace().getTypeCode()==Points) {
1001     Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1002     } else {
1003     throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1004     }
1005     break;
1006 phornby 1628 case(ContactElementsZero):
1007     case(ContactElementsOne):
1008 gross 2385 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1009     Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1010     } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1011     Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1012     } else {
1013     throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1014     }
1015     break;
1016 phornby 1628 case(ReducedContactElementsZero):
1017     case(ReducedContactElementsOne):
1018 gross 2385 if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1019     Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1020     } else {
1021     throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1022     }
1023 phornby 1628 break;
1024 gross 2385 case(DegreesOfFreedom):
1025     switch(target.getFunctionSpace().getTypeCode()) {
1026     case(ReducedDegreesOfFreedom):
1027     case(DegreesOfFreedom):
1028     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1029     break;
1030    
1031     case(Nodes):
1032     case(ReducedNodes):
1033     if (getMPISize()>1) {
1034     escript::Data temp=escript::Data(in);
1035     temp.expand();
1036     escriptDataC _in2 = temp.getDataC();
1037     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1038     } else {
1039     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1040     }
1041     break;
1042     case(Elements):
1043     case(ReducedElements):
1044     if (getMPISize()>1) {
1045     escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
1046     escriptDataC _in2 = temp.getDataC();
1047     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1048     } else {
1049     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1050     }
1051     break;
1052     case(FaceElements):
1053     case(ReducedFaceElements):
1054     if (getMPISize()>1) {
1055     escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
1056     escriptDataC _in2 = temp.getDataC();
1057     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1058    
1059     } else {
1060     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1061     }
1062     break;
1063     case(Points):
1064     if (getMPISize()>1) {
1065     escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
1066     escriptDataC _in2 = temp.getDataC();
1067     } else {
1068     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1069     }
1070     break;
1071     case(ContactElementsZero):
1072     case(ContactElementsOne):
1073     case(ReducedContactElementsZero):
1074     case(ReducedContactElementsOne):
1075     if (getMPISize()>1) {
1076     escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
1077     escriptDataC _in2 = temp.getDataC();
1078     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1079     } else {
1080     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1081     }
1082     break;
1083     default:
1084     stringstream temp;
1085     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1086     throw FinleyAdapterException(temp.str());
1087     break;
1088     }
1089     break;
1090     case(ReducedDegreesOfFreedom):
1091     switch(target.getFunctionSpace().getTypeCode()) {
1092     case(Nodes):
1093     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1094     break;
1095     case(ReducedNodes):
1096     if (getMPISize()>1) {
1097     escript::Data temp=escript::Data(in);
1098     temp.expand();
1099     escriptDataC _in2 = temp.getDataC();
1100     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1101     } else {
1102     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1103     }
1104     break;
1105     case(DegreesOfFreedom):
1106     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1107     break;
1108     case(ReducedDegreesOfFreedom):
1109     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1110     break;
1111     case(Elements):
1112     case(ReducedElements):
1113     if (getMPISize()>1) {
1114     escript::Data temp=escript::Data( in, reducedContinuousFunction(asAbstractContinuousDomain()) );
1115     escriptDataC _in2 = temp.getDataC();
1116     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1117     } else {
1118     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1119     }
1120     break;
1121     case(FaceElements):
1122     case(ReducedFaceElements):
1123     if (getMPISize()>1) {
1124     escript::Data temp=escript::Data( in, reducedContinuousFunction(asAbstractContinuousDomain()) );
1125     escriptDataC _in2 = temp.getDataC();
1126     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1127     } else {
1128     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1129     }
1130     break;
1131     case(Points):
1132     if (getMPISize()>1) {
1133     escript::Data temp=escript::Data( in, reducedContinuousFunction(asAbstractContinuousDomain()) );
1134     escriptDataC _in2 = temp.getDataC();
1135     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1136     } else {
1137     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1138     }
1139     break;
1140     case(ContactElementsZero):
1141     case(ContactElementsOne):
1142     case(ReducedContactElementsZero):
1143     case(ReducedContactElementsOne):
1144     if (getMPISize()>1) {
1145     escript::Data temp=escript::Data( in, reducedContinuousFunction(asAbstractContinuousDomain()) );
1146     escriptDataC _in2 = temp.getDataC();
1147     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1148     } else {
1149     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1150     }
1151     break;
1152     default:
1153     stringstream temp;
1154     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1155     throw FinleyAdapterException(temp.str());
1156     break;
1157     }
1158     break;
1159 phornby 1628 default:
1160     stringstream temp;
1161     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
1162     throw FinleyAdapterException(temp.str());
1163     break;
1164     }
1165     checkFinleyError();
1166 jgs 82 }
1167    
1168     //
1169     // copies the locations of sample points into x:
1170     //
1171 jgs 480 void MeshAdapter::setToX(escript::Data& arg) const
1172 jgs 82 {
1173 jfenwick 1872 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1174 phornby 1628 if (argDomain!=*this)
1175     throw FinleyAdapterException("Error - Illegal domain of data point locations");
1176     Finley_Mesh* mesh=m_finleyMesh.get();
1177     // in case of values node coordinates we can do the job directly:
1178     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1179     escriptDataC _arg=arg.getDataC();
1180     Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1181     } else {
1182     escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
1183     escriptDataC _tmp_data=tmp_data.getDataC();
1184     Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1185     // this is then interpolated onto arg:
1186     interpolateOnDomain(arg,tmp_data);
1187     }
1188     checkFinleyError();
1189 jgs 82 }
1190 jgs 149
1191 jgs 82 //
1192     // return the normal vectors at the location of data points as a Data object:
1193     //
1194 jgs 480 void MeshAdapter::setToNormal(escript::Data& normal) const
1195 jgs 82 {
1196 jfenwick 1872 /* const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1197     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1198 phornby 1628 if (normalDomain!=*this)
1199     throw FinleyAdapterException("Error - Illegal domain of normal locations");
1200     Finley_Mesh* mesh=m_finleyMesh.get();
1201     escriptDataC _normal=normal.getDataC();
1202     switch(normal.getFunctionSpace().getTypeCode()) {
1203     case(Nodes):
1204     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
1205     break;
1206     case(ReducedNodes):
1207     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
1208     break;
1209     case(Elements):
1210     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
1211     break;
1212     case(ReducedElements):
1213     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
1214     break;
1215     case (FaceElements):
1216     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1217     break;
1218     case (ReducedFaceElements):
1219     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
1220     break;
1221     case(Points):
1222     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
1223     break;
1224     case (ContactElementsOne):
1225     case (ContactElementsZero):
1226     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1227     break;
1228     case (ReducedContactElementsOne):
1229     case (ReducedContactElementsZero):
1230     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
1231     break;
1232     case(DegreesOfFreedom):
1233     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
1234     break;
1235     case(ReducedDegreesOfFreedom):
1236     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
1237     break;
1238     default:
1239 jgs 150 stringstream temp;
1240     temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
1241     throw FinleyAdapterException(temp.str());
1242 jgs 82 break;
1243 phornby 1628 }
1244     checkFinleyError();
1245 jgs 82 }
1246 jgs 149
1247 jgs 82 //
1248     // interpolates data to other domain:
1249     //
1250 jgs 480 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1251 jgs 82 {
1252 jfenwick 2087 const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1253     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1254     if (targetDomain!=this)
1255 phornby 1628 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1256 jgs 82
1257 phornby 1628 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
1258 jgs 82 }
1259 jgs 149
1260 jgs 82 //
1261     // calculates the integral of a function defined of arg:
1262     //
1263 caltinay 2778 void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
1264 jgs 82 {
1265 jfenwick 1872 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1266 phornby 1628 if (argDomain!=*this)
1267     throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1268 jgs 82
1269 phornby 1628 double blocktimer_start = blocktimer_time();
1270     Finley_Mesh* mesh=m_finleyMesh.get();
1271     escriptDataC _temp;
1272     escript::Data temp;
1273     escriptDataC _arg=arg.getDataC();
1274     switch(arg.getFunctionSpace().getTypeCode()) {
1275     case(Nodes):
1276 jfenwick 1872 temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1277 phornby 1628 _temp=temp.getDataC();
1278     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1279     break;
1280     case(ReducedNodes):
1281 jfenwick 1872 temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1282 phornby 1628 _temp=temp.getDataC();
1283     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1284     break;
1285     case(Elements):
1286     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1287     break;
1288     case(ReducedElements):
1289     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
1290     break;
1291     case(FaceElements):
1292     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1293     break;
1294     case(ReducedFaceElements):
1295     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
1296     break;
1297     case(Points):
1298     throw FinleyAdapterException("Error - Integral of data on points is not supported.");
1299     break;
1300     case(ContactElementsZero):
1301     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1302     break;
1303     case(ReducedContactElementsZero):
1304     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1305     break;
1306     case(ContactElementsOne):
1307     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1308     break;
1309     case(ReducedContactElementsOne):
1310     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1311     break;
1312     case(DegreesOfFreedom):
1313 jfenwick 1872 temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1314 phornby 1628 _temp=temp.getDataC();
1315     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1316     break;
1317     case(ReducedDegreesOfFreedom):
1318 jfenwick 1872 temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1319 phornby 1628 _temp=temp.getDataC();
1320     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1321     break;
1322     default:
1323     stringstream temp;
1324     temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1325     throw FinleyAdapterException(temp.str());
1326     break;
1327     }
1328     checkFinleyError();
1329     blocktimer_increment("integrate()", blocktimer_start);
1330 jgs 82 }
1331 jgs 149
1332 jgs 82 //
1333     // calculates the gradient of arg:
1334     //
1335 jgs 480 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1336 jgs 82 {
1337 jfenwick 1872 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1338 phornby 1628 if (argDomain!=*this)
1339     throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1340 jfenwick 1872 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1341 phornby 1628 if (gradDomain!=*this)
1342     throw FinleyAdapterException("Error - Illegal domain of gradient");
1343 jgs 82
1344 phornby 1628 Finley_Mesh* mesh=m_finleyMesh.get();
1345     escriptDataC _grad=grad.getDataC();
1346     escriptDataC nodeDataC;
1347     escript::Data temp;
1348     if (getMPISize()>1) {
1349 ksteube 1312 if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1350 phornby 1628 temp=escript::Data( arg, continuousFunction(asAbstractContinuousDomain()) );
1351     nodeDataC = temp.getDataC();
1352 ksteube 1312 } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1353 phornby 1628 temp=escript::Data( arg, reducedContinuousFunction(asAbstractContinuousDomain()) );
1354     nodeDataC = temp.getDataC();
1355 ksteube 1312 } else {
1356 phornby 1628 nodeDataC = arg.getDataC();
1357 ksteube 1312 }
1358 phornby 1628 } else {
1359     nodeDataC = arg.getDataC();
1360     }
1361     switch(grad.getFunctionSpace().getTypeCode()) {
1362     case(Nodes):
1363     throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
1364     break;
1365     case(ReducedNodes):
1366     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
1367     break;
1368     case(Elements):
1369     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1370     break;
1371     case(ReducedElements):
1372     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
1373     break;
1374     case(FaceElements):
1375     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1376     break;
1377     case(ReducedFaceElements):
1378     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
1379     break;
1380     case(Points):
1381     throw FinleyAdapterException("Error - Gradient at points is not supported.");
1382     break;
1383     case(ContactElementsZero):
1384     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1385     break;
1386     case(ReducedContactElementsZero):
1387     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1388     break;
1389     case(ContactElementsOne):
1390     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1391     break;
1392     case(ReducedContactElementsOne):
1393     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
1394     break;
1395     case(DegreesOfFreedom):
1396     throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
1397     break;
1398     case(ReducedDegreesOfFreedom):
1399     throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
1400     break;
1401     default:
1402     stringstream temp;
1403     temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1404     throw FinleyAdapterException(temp.str());
1405     break;
1406     }
1407     checkFinleyError();
1408 jgs 82 }
1409 jgs 149
1410 jgs 82 //
1411     // returns the size of elements:
1412     //
1413 jgs 480 void MeshAdapter::setToSize(escript::Data& size) const
1414 jgs 82 {
1415 phornby 1628 Finley_Mesh* mesh=m_finleyMesh.get();
1416     escriptDataC tmp=size.getDataC();
1417     switch(size.getFunctionSpace().getTypeCode()) {
1418     case(Nodes):
1419     throw FinleyAdapterException("Error - Size of nodes is not supported.");
1420     break;
1421     case(ReducedNodes):
1422     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
1423     break;
1424     case(Elements):
1425     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1426     break;
1427     case(ReducedElements):
1428     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
1429     break;
1430     case(FaceElements):
1431     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1432     break;
1433     case(ReducedFaceElements):
1434     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
1435     break;
1436     case(Points):
1437     throw FinleyAdapterException("Error - Size of point elements is not supported.");
1438     break;
1439     case(ContactElementsZero):
1440     case(ContactElementsOne):
1441     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1442     break;
1443     case(ReducedContactElementsZero):
1444     case(ReducedContactElementsOne):
1445     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
1446     break;
1447     case(DegreesOfFreedom):
1448     throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
1449     break;
1450     case(ReducedDegreesOfFreedom):
1451     throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
1452     break;
1453     default:
1454     stringstream temp;
1455     temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1456     throw FinleyAdapterException(temp.str());
1457     break;
1458     }
1459     checkFinleyError();
1460 jgs 82 }
1461 jgs 149
1462 caltinay 2151 //
1463     // sets the location of nodes
1464     //
1465 jgs 480 void MeshAdapter::setNewX(const escript::Data& new_x)
1466 jgs 82 {
1467 phornby 1628 Finley_Mesh* mesh=m_finleyMesh.get();
1468     escriptDataC tmp;
1469 jfenwick 1872 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1470 phornby 1628 if (newDomain!=*this)
1471     throw FinleyAdapterException("Error - Illegal domain of new point locations");
1472 gross 2518 if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {
1473     tmp = new_x.getDataC();
1474     Finley_Mesh_setCoordinates(mesh,&tmp);
1475     } else {
1476     escript::Data new_x_inter=escript::Data( new_x, continuousFunction(asAbstractContinuousDomain()) );
1477     tmp = new_x_inter.getDataC();
1478     Finley_Mesh_setCoordinates(mesh,&tmp);
1479     }
1480 phornby 1628 checkFinleyError();
1481 jgs 82 }
1482 jgs 149
1483 caltinay 2151 //
1484     // Helper for the save* methods. Extracts optional data variable names and the
1485     // corresponding pointers from python dictionary. Caller must free arrays.
1486     //
1487     void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1488 jgs 153 {
1489 caltinay 2151 numData = boost::python::extract<int>(arg.attr("__len__")());
1490 phornby 1628 /* win32 refactor */
1491 caltinay 2151 names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1492     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1493     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
1494 jgs 153
1495 phornby 1628 boost::python::list keys=arg.keys();
1496 caltinay 2151 for (int i=0; i<numData; ++i) {
1497 caltinay 2778 string n=boost::python::extract<string>(keys[i]);
1498 phornby 1628 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1499 jfenwick 1872 if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1500 caltinay 2151 throw FinleyAdapterException("Error: Data must be defined on same Domain");
1501     data[i] = d.getDataC();
1502     dataPtr[i] = &(data[i]);
1503 caltinay 2150 names[i] = TMPMEMALLOC(n.length()+1, char);
1504     strcpy(names[i], n.c_str());
1505 phornby 1628 }
1506 caltinay 2151 }
1507    
1508     //
1509     // saves mesh and optionally data arrays in openDX format
1510     //
1511 caltinay 2778 void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
1512 caltinay 2151 {
1513     int num_data;
1514     char **names;
1515     escriptDataC *data;
1516     escriptDataC **ptr_data;
1517    
1518     extractArgsFromDict(arg, num_data, names, data, ptr_data);
1519     Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1520 phornby 1628 checkFinleyError();
1521    
1522     /* win32 refactor */
1523     TMPMEMFREE(data);
1524     TMPMEMFREE(ptr_data);
1525 caltinay 2150 for(int i=0; i<num_data; i++)
1526 phornby 1628 {
1527     TMPMEMFREE(names[i]);
1528     }
1529     TMPMEMFREE(names);
1530 woo409 757
1531 phornby 1628 return;
1532 jgs 82 }
1533 jgs 149
1534 caltinay 2151 //
1535     // saves mesh and optionally data arrays in VTK format
1536     //
1537 caltinay 2778 void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg, const string& metadata, const string& metadata_schema) const
1538 jgs 153 {
1539 caltinay 2151 int num_data;
1540     char **names;
1541     escriptDataC *data;
1542     escriptDataC **ptr_data;
1543 jgs 153
1544 caltinay 2151 extractArgsFromDict(arg, num_data, names, data, ptr_data);
1545 gross 2421 Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1546 caltinay 2151 checkFinleyError();
1547 woo409 757
1548 phornby 1628 /* win32 refactor */
1549     TMPMEMFREE(data);
1550     TMPMEMFREE(ptr_data);
1551 caltinay 2151 for(int i=0; i<num_data; i++)
1552 phornby 1628 {
1553     TMPMEMFREE(names[i]);
1554     }
1555     TMPMEMFREE(names);
1556 caltinay 2151 }
1557 woo409 757
1558 jfenwick 2642 bool MeshAdapter::ownSample(int fs_code, index_t id) const
1559     {
1560 jfenwick 2644 #ifdef PASO_MPI
1561 jfenwick 2642 index_t myFirstNode=0, myLastNode=0, k=0;
1562     index_t* globalNodeIndex=0;
1563     Finley_Mesh* mesh_p=m_finleyMesh.get();
1564     if (fs_code == FINLEY_REDUCED_NODES)
1565     {
1566     myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1567     myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1568     globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1569     }
1570     else
1571     {
1572     myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1573     myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1574     globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1575     }
1576     k=globalNodeIndex[id];
1577     return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1578 jfenwick 2644 #endif
1579     return true;
1580 jfenwick 2642 }
1581    
1582    
1583    
1584 caltinay 2151 //
1585     // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1586     //
1587 jgs 82 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1588 phornby 1628 const int row_blocksize,
1589     const escript::FunctionSpace& row_functionspace,
1590     const int column_blocksize,
1591     const escript::FunctionSpace& column_functionspace,
1592     const int type) const
1593 jgs 82 {
1594 phornby 1628 int reduceRowOrder=0;
1595     int reduceColOrder=0;
1596     // is the domain right?
1597 jfenwick 1872 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1598 phornby 1628 if (row_domain!=*this)
1599     throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1600 jfenwick 1872 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1601 phornby 1628 if (col_domain!=*this)
1602     throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1603     // is the function space type right
1604     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1605     reduceRowOrder=0;
1606     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1607     reduceRowOrder=1;
1608     } else {
1609     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1610     }
1611     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1612     reduceColOrder=0;
1613     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1614     reduceColOrder=1;
1615     } else {
1616     throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1617     }
1618     // generate matrix:
1619    
1620     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1621     checkFinleyError();
1622     Paso_SystemMatrix* fsystemMatrix;
1623     int trilinos = 0;
1624     if (trilinos) {
1625 ksteube 1312 #ifdef TRILINOS
1626     /* Allocation Epetra_VrbMatrix here */
1627     #endif
1628 phornby 1628 }
1629     else {
1630 gross 2551 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1631 phornby 1628 }
1632     checkPasoError();
1633     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1634     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1635 jgs 82 }
1636 caltinay 2151
1637     //
1638 gross 1366 // creates a TransportProblemAdapter
1639 caltinay 2151 //
1640 gross 1366 TransportProblemAdapter MeshAdapter::newTransportProblem(
1641 phornby 1628 const double theta,
1642     const int blocksize,
1643     const escript::FunctionSpace& functionspace,
1644     const int type) const
1645 gross 1366 {
1646 phornby 1628 int reduceOrder=0;
1647     // is the domain right?
1648 jfenwick 1872 const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1649 phornby 1628 if (domain!=*this)
1650     throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1651     // is the function space type right
1652     if (functionspace.getTypeCode()==DegreesOfFreedom) {
1653     reduceOrder=0;
1654     } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1655     reduceOrder=1;
1656     } else {
1657     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1658     }
1659     // generate matrix:
1660    
1661     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1662     checkFinleyError();
1663     Paso_FCTransportProblem* transportProblem;
1664     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);
1665     checkPasoError();
1666     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1667     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);
1668 gross 1366 }
1669 jgs 149
1670 jgs 82 //
1671     // vtkObject MeshAdapter::createVtkObject() const
1672     // TODO:
1673     //
1674 jgs 149
1675 jgs 82 //
1676     // returns true if data at the atom_type is considered as being cell centered:
1677     bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1678     {
1679 phornby 1628 switch(functionSpaceCode) {
1680     case(Nodes):
1681     case(DegreesOfFreedom):
1682     case(ReducedDegreesOfFreedom):
1683     return false;
1684     break;
1685     case(Elements):
1686     case(FaceElements):
1687     case(Points):
1688     case(ContactElementsZero):
1689     case(ContactElementsOne):
1690     case(ReducedElements):
1691     case(ReducedFaceElements):
1692     case(ReducedContactElementsZero):
1693     case(ReducedContactElementsOne):
1694     return true;
1695     break;
1696     default:
1697     stringstream temp;
1698     temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1699     throw FinleyAdapterException(temp.str());
1700     break;
1701     }
1702     checkFinleyError();
1703     return false;
1704 jgs 82 }
1705 jgs 149
1706 jfenwick 2635 bool
1707 caltinay 2778 MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1708 jfenwick 2635 {
1709     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1710     class 1: DOF <-> Nodes
1711     class 2: ReducedDOF <-> ReducedNodes
1712     class 3: Points
1713     class 4: Elements
1714     class 5: ReducedElements
1715     class 6: FaceElements
1716     class 7: ReducedFaceElements
1717     class 8: ContactElementZero <-> ContactElementOne
1718     class 9: ReducedContactElementZero <-> ReducedContactElementOne
1719    
1720     There is also a set of lines. Interpolation is possible down a line but not between lines.
1721     class 1 and 2 belong to all lines so aren't considered.
1722     line 0: class 3
1723     line 1: class 4,5
1724     line 2: class 6,7
1725     line 3: class 8,9
1726    
1727     For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1728     eg hasnodes is true if we have at least one instance of Nodes.
1729     */
1730     if (fs.size()==0)
1731     {
1732 caltinay 2778 return false;
1733 jfenwick 2635 }
1734 caltinay 2778 vector<int> hasclass(10);
1735     vector<int> hasline(4);
1736 jfenwick 2635 bool hasnodes=false;
1737     bool hasrednodes=false;
1738     bool hascez=false;
1739     bool hasrcez=false;
1740     for (int i=0;i<fs.size();++i)
1741     {
1742     switch(fs[i])
1743     {
1744     case(Nodes): hasnodes=true; // no break is deliberate
1745     case(DegreesOfFreedom):
1746     hasclass[1]=1;
1747     break;
1748     case(ReducedNodes): hasrednodes=true; // no break is deliberate
1749     case(ReducedDegreesOfFreedom):
1750     hasclass[2]=1;
1751     break;
1752     case(Points):
1753     hasline[0]=1;
1754     hasclass[3]=1;
1755     break;
1756     case(Elements):
1757     hasclass[4]=1;
1758     hasline[1]=1;
1759     break;
1760     case(ReducedElements):
1761     hasclass[5]=1;
1762     hasline[1]=1;
1763     break;
1764     case(FaceElements):
1765     hasclass[6]=1;
1766     hasline[2]=1;
1767     break;
1768     case(ReducedFaceElements):
1769     hasclass[7]=1;
1770     hasline[2]=1;
1771     break;
1772     case(ContactElementsZero): hascez=true; // no break is deliberate
1773     case(ContactElementsOne):
1774     hasclass[8]=1;
1775     hasline[3]=1;
1776     break;
1777     case(ReducedContactElementsZero): hasrcez=true; // no break is deliberate
1778     case(ReducedContactElementsOne):
1779     hasclass[9]=1;
1780     hasline[3]=1;
1781     break;
1782     default:
1783     return false;
1784     }
1785     }
1786     int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1787     // fail if we have more than one leaf group
1788    
1789     if (totlines>1)
1790     {
1791     return false; // there are at least two branches we can't interpolate between
1792     }
1793     else if (totlines==1)
1794     {
1795     if (hasline[0]==1) // we have points
1796     {
1797     resultcode=Points;
1798     }
1799     else if (hasline[1]==1)
1800     {
1801     if (hasclass[5]==1)
1802     {
1803     resultcode=ReducedElements;
1804     }
1805     else
1806     {
1807     resultcode=Elements;
1808     }
1809     }
1810     else if (hasline[2]==1)
1811     {
1812 jfenwick 2644 if (hasclass[7]==1)
1813 jfenwick 2635 {
1814     resultcode=ReducedFaceElements;
1815     }
1816     else
1817     {
1818     resultcode=FaceElements;
1819     }
1820     }
1821     else // so we must be in line3
1822     {
1823     if (hasclass[9]==1)
1824     {
1825     // need something from class 9
1826     resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1827     }
1828     else
1829     {
1830     // something from class 8
1831     resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1832     }
1833     }
1834     }
1835     else // totlines==0
1836     {
1837     if (hasclass[2]==1)
1838     {
1839     // something from class 2
1840     resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1841     }
1842     else
1843     { // something from class 1
1844     resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1845     }
1846     }
1847     return true;
1848     }
1849    
1850 jgs 82 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1851     {
1852 phornby 1628 switch(functionSpaceType_source) {
1853     case(Nodes):
1854 jfenwick 2635 switch(functionSpaceType_target) {
1855     case(Nodes):
1856     case(ReducedNodes):
1857     case(ReducedDegreesOfFreedom):
1858     case(DegreesOfFreedom):
1859     case(Elements):
1860     case(ReducedElements):
1861     case(FaceElements):
1862     case(ReducedFaceElements):
1863     case(Points):
1864     case(ContactElementsZero):
1865     case(ReducedContactElementsZero):
1866     case(ContactElementsOne):
1867     case(ReducedContactElementsOne):
1868     return true;
1869     default:
1870     stringstream temp;
1871     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1872     throw FinleyAdapterException(temp.str());
1873 phornby 1628 }
1874     break;
1875     case(ReducedNodes):
1876 jfenwick 2635 switch(functionSpaceType_target) {
1877     case(ReducedNodes):
1878     case(ReducedDegreesOfFreedom):
1879     case(Elements):
1880     case(ReducedElements):
1881     case(FaceElements):
1882     case(ReducedFaceElements):
1883     case(Points):
1884     case(ContactElementsZero):
1885     case(ReducedContactElementsZero):
1886     case(ContactElementsOne):
1887     case(ReducedContactElementsOne):
1888     return true;
1889     case(Nodes):
1890     case(DegreesOfFreedom):
1891     return false;
1892     default:
1893     stringstream temp;
1894     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1895     throw FinleyAdapterException(temp.str());
1896 phornby 1628 }
1897     break;
1898     case(Elements):
1899 jfenwick 2635 if (functionSpaceType_target==Elements) {
1900     return true;
1901     } else if (functionSpaceType_target==ReducedElements) {
1902     return true;
1903     } else {
1904     return false;
1905     }
1906 phornby 1628 case(ReducedElements):
1907 jfenwick 2635 if (functionSpaceType_target==ReducedElements) {
1908     return true;
1909     } else {
1910     return false;
1911     }
1912 phornby 1628 case(FaceElements):
1913 jfenwick 2635 if (functionSpaceType_target==FaceElements) {
1914     return true;
1915     } else if (functionSpaceType_target==ReducedFaceElements) {
1916     return true;
1917     } else {
1918     return false;
1919     }
1920 phornby 1628 case(ReducedFaceElements):
1921 jfenwick 2635 if (functionSpaceType_target==ReducedFaceElements) {
1922     return true;
1923     } else {
1924     return false;
1925     }
1926 phornby 1628 case(Points):
1927 jfenwick 2635 if (functionSpaceType_target==Points) {
1928     return true;
1929     } else {
1930     return false;
1931     }
1932 phornby 1628 case(ContactElementsZero):
1933     case(ContactElementsOne):
1934 jfenwick 2635 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1935     return true;
1936     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1937     return true;
1938     } else {
1939     return false;
1940     }
1941 phornby 1628 case(ReducedContactElementsZero):
1942     case(ReducedContactElementsOne):
1943 jfenwick 2635 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1944     return true;
1945     } else {
1946     return false;
1947     }
1948 phornby 1628 case(DegreesOfFreedom):
1949 jfenwick 2635 switch(functionSpaceType_target) {
1950     case(ReducedDegreesOfFreedom):
1951     case(DegreesOfFreedom):
1952     case(Nodes):
1953     case(ReducedNodes):
1954     case(Elements):
1955     case(ReducedElements):
1956     case(Points):
1957     case(FaceElements):
1958     case(ReducedFaceElements):
1959     case(ContactElementsZero):
1960     case(ReducedContactElementsZero):
1961     case(ContactElementsOne):
1962     case(ReducedContactElementsOne):
1963     return true;
1964     default:
1965     stringstream temp;
1966     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1967     throw FinleyAdapterException(temp.str());
1968     }
1969     break;
1970 phornby 1628 case(ReducedDegreesOfFreedom):
1971     switch(functionSpaceType_target) {
1972 jfenwick 2635 case(ReducedDegreesOfFreedom):
1973     case(ReducedNodes):
1974     case(Elements):
1975     case(ReducedElements):
1976     case(FaceElements):
1977     case(ReducedFaceElements):
1978     case(Points):
1979     case(ContactElementsZero):
1980     case(ReducedContactElementsZero):
1981     case(ContactElementsOne):
1982     case(ReducedContactElementsOne):
1983     return true;
1984     case(Nodes):
1985     case(DegreesOfFreedom):
1986     return false;
1987     default:
1988     stringstream temp;
1989     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1990     throw FinleyAdapterException(temp.str());
1991     }
1992     break;
1993 phornby 1628 default:
1994     stringstream temp;
1995     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1996     throw FinleyAdapterException(temp.str());
1997     break;
1998     }
1999     checkFinleyError();
2000     return false;
2001 jgs 82 }
2002 jgs 149
2003 jgs 82 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
2004     {
2005 phornby 1628 return false;
2006 jgs 82 }
2007 jgs 104
2008 jgs 121 bool MeshAdapter::operator==(const AbstractDomain& other) const
2009 jgs 82 {
2010 phornby 1628 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
2011     if (temp!=0) {
2012     return (m_finleyMesh==temp->m_finleyMesh);
2013     } else {
2014     return false;
2015     }
2016 jgs 82 }
2017    
2018 jgs 121 bool MeshAdapter::operator!=(const AbstractDomain& other) const
2019 jgs 104 {
2020 phornby 1628 return !(operator==(other));
2021 jgs 104 }
2022    
2023 gross 1859 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2024 jgs 102 {
2025 gross 2315 Finley_Mesh* mesh=m_finleyMesh.get();
2026     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2027 jgs 150 checkPasoError();
2028 jgs 102 return out;
2029     }
2030 gross 1859 int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2031     {
2032 gross 2315 Finley_Mesh* mesh=m_finleyMesh.get();
2033     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2034 gross 1859 checkPasoError();
2035     return out;
2036     }
2037 jgs 149
2038 jgs 480 escript::Data MeshAdapter::getX() const
2039 jgs 102 {
2040 phornby 1628 return continuousFunction(asAbstractContinuousDomain()).getX();
2041 jgs 102 }
2042 jgs 149
2043 jgs 480 escript::Data MeshAdapter::getNormal() const
2044 jgs 102 {
2045 phornby 1628 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
2046 jgs 102 }
2047 jgs 149
2048 jgs 480 escript::Data MeshAdapter::getSize() const
2049 jgs 102 {
2050 jfenwick 1872 return escript::function(asAbstractContinuousDomain()).getSize();
2051 jgs 102 }
2052    
2053 jfenwick 2487 const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2054 gross 1062 {
2055 phornby 1628 int *out = NULL;
2056     Finley_Mesh* mesh=m_finleyMesh.get();
2057     switch (functionSpaceType) {
2058     case(Nodes):
2059     out=mesh->Nodes->Id;
2060     break;
2061     case(ReducedNodes):
2062     out=mesh->Nodes->reducedNodesId;
2063     break;
2064     case(Elements):
2065     out=mesh->Elements->Id;
2066     break;
2067     case(ReducedElements):
2068     out=mesh->Elements->Id;
2069     break;
2070     case(FaceElements):
2071     out=mesh->FaceElements->Id;
2072     break;
2073     case(ReducedFaceElements):
2074     out=mesh->FaceElements->Id;
2075     break;
2076     case(Points):
2077     out=mesh->Points->Id;
2078     break;
2079     case(ContactElementsZero):
2080     out=mesh->ContactElements->Id;
2081     break;
2082     case(ReducedContactElementsZero):
2083     out=mesh->ContactElements->Id;
2084     break;
2085     case(ContactElementsOne):
2086     out=mesh->ContactElements->Id;
2087     break;
2088     case(ReducedContactElementsOne):
2089     out=mesh->ContactElements->Id;
2090     break;
2091     case(DegreesOfFreedom):
2092     out=mesh->Nodes->degreesOfFreedomId;
2093     break;
2094     case(ReducedDegreesOfFreedom):
2095     out=mesh->Nodes->reducedDegreesOfFreedomId;
2096     break;
2097     default:
2098 gross 1062 stringstream temp;
2099     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
2100     throw FinleyAdapterException(temp.str());
2101     break;
2102 phornby 1628 }
2103     return out;
2104 gross 1062 }
2105 jgs 110 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
2106     {
2107 phornby 1628 int out=0;
2108     Finley_Mesh* mesh=m_finleyMesh.get();
2109     switch (functionSpaceType) {
2110     case(Nodes):
2111     out=mesh->Nodes->Tag[sampleNo];
2112     break;
2113     case(ReducedNodes):
2114     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
2115     break;
2116     case(Elements):
2117     out=mesh->Elements->Tag[sampleNo];
2118     break;
2119     case(ReducedElements):
2120     out=mesh->Elements->Tag[sampleNo];
2121     break;
2122     case(FaceElements):
2123     out=mesh->FaceElements->Tag[sampleNo];
2124     break;
2125     case(ReducedFaceElements):
2126     out=mesh->FaceElements->Tag[sampleNo];
2127     break;
2128     case(Points):
2129     out=mesh->Points->Tag[sampleNo];
2130     break;
2131     case(ContactElementsZero):
2132     out=mesh->ContactElements->Tag[sampleNo];
2133     break;
2134     case(ReducedContactElementsZero):
2135     out=mesh->ContactElements->Tag[sampleNo];
2136     break;
2137     case(ContactElementsOne):
2138     out=mesh->ContactElements->Tag[sampleNo];
2139     break;
2140     case(ReducedContactElementsOne):
2141     out=mesh->ContactElements->Tag[sampleNo];
2142     break;
2143     case(DegreesOfFreedom):
2144     throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
2145     break;
2146     case(ReducedDegreesOfFreedom):
2147     throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
2148     break;
2149     default:
2150 &nb