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