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