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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (hide annotations)
Mon Oct 11 01:48:14 2010 UTC (8 years, 9 months ago) by jfenwick
File size: 85719 byte(s)
Merging dudley and scons updates from branches

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