/[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 1748 - (hide annotations)
Wed Sep 3 06:10:39 2008 UTC (10 years, 9 months ago) by ksteube
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 84232 byte(s)
MPI parallelism for Data().dump and load.  Use multiple NetCDF
files, one file per MPI process

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26