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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2149 - (hide annotations)
Wed Dec 10 05:08:17 2008 UTC (10 years, 4 months ago) by caltinay
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 85399 byte(s)
Fixed a typo in error message.

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