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

Contents of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5148 - (show annotations)
Mon Sep 15 01:25:23 2014 UTC (4 years, 11 months ago) by caltinay
File size: 91360 byte(s)
Merging ripley diagonal storage + CUDA support into trunk.
Options file version has been incremented due to new options
'cuda' and 'nvccflags'.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26