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

Contents of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5122 - (show annotations)
Thu Aug 21 04:32:39 2014 UTC (4 years, 6 months ago) by caltinay
File size: 91592 byte(s)
Fast-forward to latest trunk to be able to use both Paso or cusp.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26