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