/[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 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (11 years, 6 months ago) by ksteube
File size: 58752 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

1
2 /* $Id$ */
3
4 /*******************************************************
5 *
6 * Copyright 2003-2007 by ACceSS MNRF
7 * Copyright 2007 by University of Queensland
8 *
9 * http://esscc.uq.edu.au
10 * Primary Business: Queensland, Australia
11 * Licensed under the Open Software License version 3.0
12 * http://www.opensource.org/licenses/osl-3.0.php
13 *
14 *******************************************************/
15
16 #include "MeshAdapter.h"
17 #include "escript/Data.h"
18 #include "escript/DataFactory.h"
19 extern "C" {
20 #include "escript/blocktimer.h"
21 }
22 #include <vector>
23
24 using namespace std;
25 using namespace escript;
26
27 namespace finley {
28
29 //
30 // define the static constants
31 MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
32 const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
33 const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
34 const int MeshAdapter::Nodes=FINLEY_NODES;
35 const int MeshAdapter::ReducedNodes=FINLEY_REDUCED_NODES;
36 const int MeshAdapter::Elements=FINLEY_ELEMENTS;
37 const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;
38 const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
39 const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;
40 const int MeshAdapter::Points=FINLEY_POINTS;
41 const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
42 const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;
43 const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
44 const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
45
46 MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
47 {
48 setFunctionSpaceTypeNames();
49 //
50 // need to use a null_deleter as Finley_Mesh_free deletes the pointer
51 // for us.
52 m_finleyMesh.reset(finleyMesh,null_deleter());
53 }
54
55 //
56 // The copy constructor should just increment the use count
57 MeshAdapter::MeshAdapter(const MeshAdapter& in):
58 m_finleyMesh(in.m_finleyMesh)
59 {
60 setFunctionSpaceTypeNames();
61 }
62
63 MeshAdapter::~MeshAdapter()
64 {
65 //
66 // I hope the case for the pointer being zero has been taken care of.
67 // cout << "In MeshAdapter destructor." << endl;
68 if (m_finleyMesh.unique()) {
69 Finley_Mesh_free(m_finleyMesh.get());
70 }
71 }
72
73 int MeshAdapter::getMPISize() const
74 {
75 return m_finleyMesh.get()->MPIInfo->size;
76 }
77 int MeshAdapter::getMPIRank() const
78 {
79 return m_finleyMesh.get()->MPIInfo->rank;
80 }
81
82
83 Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
84 return m_finleyMesh.get();
85 }
86
87 void MeshAdapter::write(const std::string& fileName) const
88 {
89 char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
90 strcpy(fName,fileName.c_str());
91 Finley_Mesh_write(m_finleyMesh.get(),fName);
92 checkFinleyError();
93 TMPMEMFREE(fName);
94 }
95
96 void MeshAdapter::dump(const std::string& fileName) const
97 {
98 char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
99 strcpy(fName,fileName.c_str());
100 Finley_Mesh_dump(m_finleyMesh.get(),fName);
101 checkFinleyError();
102 TMPMEMFREE(fName);
103 }
104
105 string MeshAdapter::getDescription() const
106 {
107 return "FinleyMesh";
108 }
109
110 string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
111 {
112 FunctionSpaceNamesMapType::iterator loc;
113 loc=m_functionSpaceTypeNames.find(functionSpaceType);
114 if (loc==m_functionSpaceTypeNames.end()) {
115 return "Invalid function space type code.";
116 } else {
117 return loc->second;
118 }
119 }
120
121 bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
122 {
123 FunctionSpaceNamesMapType::iterator loc;
124 loc=m_functionSpaceTypeNames.find(functionSpaceType);
125 return (loc!=m_functionSpaceTypeNames.end());
126 }
127
128 void MeshAdapter::setFunctionSpaceTypeNames()
129 {
130 m_functionSpaceTypeNames.insert
131 (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));
132 m_functionSpaceTypeNames.insert
133 (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));
134 m_functionSpaceTypeNames.insert
135 (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
136 m_functionSpaceTypeNames.insert
137 (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
138 m_functionSpaceTypeNames.insert
139 (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
140 m_functionSpaceTypeNames.insert
141 (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
142 m_functionSpaceTypeNames.insert
143 (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
144 m_functionSpaceTypeNames.insert
145 (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
146 m_functionSpaceTypeNames.insert
147 (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
148 m_functionSpaceTypeNames.insert
149 (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
150 m_functionSpaceTypeNames.insert
151 (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
152 m_functionSpaceTypeNames.insert
153 (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
154 m_functionSpaceTypeNames.insert
155 (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
156 }
157
158 int MeshAdapter::getContinuousFunctionCode() const
159 {
160 return Nodes;
161 }
162 int MeshAdapter::getReducedContinuousFunctionCode() const
163 {
164 return ReducedNodes;
165 }
166
167 int MeshAdapter::getFunctionCode() const
168 {
169 return Elements;
170 }
171 int MeshAdapter::getReducedFunctionCode() const
172 {
173 return ReducedElements;
174 }
175
176 int MeshAdapter::getFunctionOnBoundaryCode() const
177 {
178 return FaceElements;
179 }
180 int MeshAdapter::getReducedFunctionOnBoundaryCode() const
181 {
182 return ReducedFaceElements;
183 }
184
185 int MeshAdapter::getFunctionOnContactZeroCode() const
186 {
187 return ContactElementsZero;
188 }
189 int MeshAdapter::getReducedFunctionOnContactZeroCode() const
190 {
191 return ReducedContactElementsZero;
192 }
193
194 int MeshAdapter::getFunctionOnContactOneCode() const
195 {
196 return ContactElementsOne;
197 }
198 int MeshAdapter::getReducedFunctionOnContactOneCode() const
199 {
200 return ReducedContactElementsOne;
201 }
202
203 int MeshAdapter::getSolutionCode() const
204 {
205 return DegreesOfFreedom;
206 }
207
208 int MeshAdapter::getReducedSolutionCode() const
209 {
210 return ReducedDegreesOfFreedom;
211 }
212
213 int MeshAdapter::getDiracDeltaFunctionCode() const
214 {
215 return Points;
216 }
217
218 //
219 // return the spatial dimension of the Mesh:
220 //
221 int MeshAdapter::getDim() const
222 {
223 int numDim=Finley_Mesh_getDim(m_finleyMesh.get());
224 checkFinleyError();
225 return numDim;
226 }
227
228 //
229 // return the number of data points per sample and the number of samples
230 // needed to represent data on a parts of the mesh.
231 //
232 pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
233 {
234 int numDataPointsPerSample=0;
235 int numSamples=0;
236 Finley_Mesh* mesh=m_finleyMesh.get();
237 switch (functionSpaceCode) {
238 case(Nodes):
239 numDataPointsPerSample=1;
240 numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
241 break;
242 case(ReducedNodes):
243 numDataPointsPerSample=1;
244 numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
245 break;
246 case(Elements):
247 if (mesh->Elements!=NULL) {
248 numSamples=mesh->Elements->numElements;
249 numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
250 }
251 break;
252 case(ReducedElements):
253 if (mesh->Elements!=NULL) {
254 numSamples=mesh->Elements->numElements;
255 numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;
256 }
257 break;
258 case(FaceElements):
259 if (mesh->FaceElements!=NULL) {
260 numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
261 numSamples=mesh->FaceElements->numElements;
262 }
263 break;
264 case(ReducedFaceElements):
265 if (mesh->FaceElements!=NULL) {
266 numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;
267 numSamples=mesh->FaceElements->numElements;
268 }
269 break;
270 case(Points):
271 if (mesh->Points!=NULL) {
272 numDataPointsPerSample=1;
273 numSamples=mesh->Points->numElements;
274 }
275 break;
276 case(ContactElementsZero):
277 if (mesh->ContactElements!=NULL) {
278 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
279 numSamples=mesh->ContactElements->numElements;
280 }
281 break;
282 case(ReducedContactElementsZero):
283 if (mesh->ContactElements!=NULL) {
284 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
285 numSamples=mesh->ContactElements->numElements;
286 }
287 break;
288 case(ContactElementsOne):
289 if (mesh->ContactElements!=NULL) {
290 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
291 numSamples=mesh->ContactElements->numElements;
292 }
293 break;
294 case(ReducedContactElementsOne):
295 if (mesh->ContactElements!=NULL) {
296 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
297 numSamples=mesh->ContactElements->numElements;
298 }
299 break;
300 case(DegreesOfFreedom):
301 if (mesh->Nodes!=NULL) {
302 numDataPointsPerSample=1;
303 numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
304 }
305 break;
306 case(ReducedDegreesOfFreedom):
307 if (mesh->Nodes!=NULL) {
308 numDataPointsPerSample=1;
309 numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
310 }
311 break;
312 default:
313 stringstream temp;
314 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
315 throw FinleyAdapterException(temp.str());
316 break;
317 }
318 return pair<int,int>(numDataPointsPerSample,numSamples);
319 }
320
321 //
322 // adds linear PDE of second order into a given stiffness matrix and right hand side:
323 //
324 void MeshAdapter::addPDEToSystem(
325 SystemMatrixAdapter& mat, escript::Data& rhs,
326 const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,const escript::Data& X,const escript::Data& Y,
327 const escript::Data& d, const escript::Data& y,
328 const escript::Data& d_contact,const escript::Data& y_contact) const
329 {
330 escriptDataC _rhs=rhs.getDataC();
331 escriptDataC _A =A.getDataC();
332 escriptDataC _B=B.getDataC();
333 escriptDataC _C=C.getDataC();
334 escriptDataC _D=D.getDataC();
335 escriptDataC _X=X.getDataC();
336 escriptDataC _Y=Y.getDataC();
337 escriptDataC _d=d.getDataC();
338 escriptDataC _y=y.getDataC();
339 escriptDataC _d_contact=d_contact.getDataC();
340 escriptDataC _y_contact=y_contact.getDataC();
341
342 Finley_Mesh* mesh=m_finleyMesh.get();
343
344 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
345 checkFinleyError();
346
347 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
348 checkFinleyError();
349
350 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
351 checkFinleyError();
352 }
353
354 void MeshAdapter::addPDEToLumpedSystem(
355 escript::Data& mat,
356 const escript::Data& D,
357 const escript::Data& d) const
358 {
359 escriptDataC _mat=mat.getDataC();
360 escriptDataC _D=D.getDataC();
361 escriptDataC _d=d.getDataC();
362
363 Finley_Mesh* mesh=m_finleyMesh.get();
364
365 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
366 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
367
368 checkFinleyError();
369 }
370
371
372 //
373 // adds linear PDE of second order into the right hand side only
374 //
375 void MeshAdapter::addPDEToRHS( escript::Data& rhs, const escript::Data& X,const escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
376 {
377 Finley_Mesh* mesh=m_finleyMesh.get();
378
379 escriptDataC _rhs=rhs.getDataC();
380 escriptDataC _X=X.getDataC();
381 escriptDataC _Y=Y.getDataC();
382 escriptDataC _y=y.getDataC();
383 escriptDataC _y_contact=y_contact.getDataC();
384
385 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
386 checkFinleyError();
387
388 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
389 checkFinleyError();
390
391 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
392 checkFinleyError();
393 }
394
395 //
396 // interpolates data between different function spaces:
397 //
398 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
399 {
400 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
401 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
402 if (inDomain!=*this)
403 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
404 if (targetDomain!=*this)
405 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
406
407 Finley_Mesh* mesh=m_finleyMesh.get();
408 escriptDataC _target=target.getDataC();
409 escriptDataC _in=in.getDataC();
410 switch(in.getFunctionSpace().getTypeCode()) {
411 case(Nodes):
412 switch(target.getFunctionSpace().getTypeCode()) {
413 case(Nodes):
414 case(ReducedNodes):
415 case(DegreesOfFreedom):
416 case(ReducedDegreesOfFreedom):
417 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
418 break;
419 case(Elements):
420 case(ReducedElements):
421 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
422 break;
423 case(FaceElements):
424 case(ReducedFaceElements):
425 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
426 break;
427 case(Points):
428 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
429 break;
430 case(ContactElementsZero):
431 case(ReducedContactElementsZero):
432 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
433 break;
434 case(ContactElementsOne):
435 case(ReducedContactElementsOne):
436 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
437 break;
438 default:
439 stringstream temp;
440 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
441 throw FinleyAdapterException(temp.str());
442 break;
443 }
444 break;
445 case(ReducedNodes):
446 switch(target.getFunctionSpace().getTypeCode()) {
447 case(Nodes):
448 case(ReducedNodes):
449 case(DegreesOfFreedom):
450 case(ReducedDegreesOfFreedom):
451 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
452 break;
453 case(Elements):
454 case(ReducedElements):
455 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
456 break;
457 case(FaceElements):
458 case(ReducedFaceElements):
459 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
460 break;
461 case(Points):
462 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
463 break;
464 case(ContactElementsZero):
465 case(ReducedContactElementsZero):
466 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
467 break;
468 case(ContactElementsOne):
469 case(ReducedContactElementsOne):
470 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
471 break;
472 default:
473 stringstream temp;
474 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
475 throw FinleyAdapterException(temp.str());
476 break;
477 }
478 break;
479 case(Elements):
480 if (target.getFunctionSpace().getTypeCode()==Elements) {
481 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
482 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
483 Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
484 } else {
485 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
486 }
487 break;
488 case(ReducedElements):
489 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
490 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
491 } else {
492 throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
493 }
494 break;
495 case(FaceElements):
496 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
497 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
498 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
499 Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
500 } else {
501 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
502 }
503 break;
504 case(ReducedFaceElements):
505 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
506 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
507 } else {
508 throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
509 }
510 break;
511 case(Points):
512 if (target.getFunctionSpace().getTypeCode()==Points) {
513 Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
514 } else {
515 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
516 }
517 break;
518 case(ContactElementsZero):
519 case(ContactElementsOne):
520 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
521 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
522 } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
523 Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
524 } else {
525 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
526 break;
527 }
528 break;
529 case(ReducedContactElementsZero):
530 case(ReducedContactElementsOne):
531 if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
532 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
533 } else {
534 throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
535 break;
536 }
537 break;
538 case(DegreesOfFreedom):
539 switch(target.getFunctionSpace().getTypeCode()) {
540 case(ReducedDegreesOfFreedom):
541 case(DegreesOfFreedom):
542 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
543 break;
544
545 case(Nodes):
546 case(ReducedNodes):
547 if (getMPISize()>1) {
548 escript::Data temp=escript::Data(in);
549 temp.expand();
550 escriptDataC _in2 = temp.getDataC();
551 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
552 } else {
553 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
554 }
555 break;
556 case(Elements):
557 case(ReducedElements):
558 if (getMPISize()>1) {
559 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
560 escriptDataC _in2 = temp.getDataC();
561 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
562 } else {
563 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
564 }
565 break;
566 case(FaceElements):
567 case(ReducedFaceElements):
568 if (getMPISize()>1) {
569 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
570 escriptDataC _in2 = temp.getDataC();
571 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
572
573 } else {
574 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
575 }
576 break;
577 case(Points):
578 if (getMPISize()>1) {
579 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
580 escriptDataC _in2 = temp.getDataC();
581 } else {
582 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
583 }
584 break;
585 case(ContactElementsZero):
586 case(ContactElementsOne):
587 case(ReducedContactElementsZero):
588 case(ReducedContactElementsOne):
589 if (getMPISize()>1) {
590 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
591 escriptDataC _in2 = temp.getDataC();
592 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
593 } else {
594 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
595 }
596 break;
597 default:
598 stringstream temp;
599 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
600 throw FinleyAdapterException(temp.str());
601 break;
602 }
603 break;
604 case(ReducedDegreesOfFreedom):
605 switch(target.getFunctionSpace().getTypeCode()) {
606 case(Nodes):
607 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
608 break;
609 case(ReducedNodes):
610 if (getMPISize()>1) {
611 escript::Data temp=escript::Data(in);
612 temp.expand();
613 escriptDataC _in2 = temp.getDataC();
614 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
615 } else {
616 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
617 }
618 break;
619 case(DegreesOfFreedom):
620 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
621 break;
622 case(ReducedDegreesOfFreedom):
623 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
624 break;
625 case(Elements):
626 case(ReducedElements):
627 if (getMPISize()>1) {
628 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
629 escriptDataC _in2 = temp.getDataC();
630 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
631 } else {
632 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
633 }
634 break;
635 case(FaceElements):
636 case(ReducedFaceElements):
637 if (getMPISize()>1) {
638 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
639 escriptDataC _in2 = temp.getDataC();
640 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
641 } else {
642 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
643 }
644 break;
645 case(Points):
646 if (getMPISize()>1) {
647 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
648 escriptDataC _in2 = temp.getDataC();
649 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
650 } else {
651 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
652 }
653 break;
654 case(ContactElementsZero):
655 case(ContactElementsOne):
656 case(ReducedContactElementsZero):
657 case(ReducedContactElementsOne):
658 if (getMPISize()>1) {
659 escript::Data temp=escript::Data( in, continuousFunction(asAbstractContinuousDomain()) );
660 escriptDataC _in2 = temp.getDataC();
661 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
662 } else {
663 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
664 }
665 break;
666 default:
667 stringstream temp;
668 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
669 throw FinleyAdapterException(temp.str());
670 break;
671 }
672 break;
673 default:
674 stringstream temp;
675 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
676 throw FinleyAdapterException(temp.str());
677 break;
678 }
679 checkFinleyError();
680 }
681
682 //
683 // copies the locations of sample points into x:
684 //
685 void MeshAdapter::setToX(escript::Data& arg) const
686 {
687 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
688 if (argDomain!=*this)
689 throw FinleyAdapterException("Error - Illegal domain of data point locations");
690 Finley_Mesh* mesh=m_finleyMesh.get();
691 // in case of values node coordinates we can do the job directly:
692 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
693 escriptDataC _arg=arg.getDataC();
694 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
695 } else {
696 escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
697 escriptDataC _tmp_data=tmp_data.getDataC();
698 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
699 // this is then interpolated onto arg:
700 interpolateOnDomain(arg,tmp_data);
701 }
702 checkFinleyError();
703 }
704
705 //
706 // return the normal vectors at the location of data points as a Data object:
707 //
708 void MeshAdapter::setToNormal(escript::Data& normal) const
709 {
710 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
711 if (normalDomain!=*this)
712 throw FinleyAdapterException("Error - Illegal domain of normal locations");
713 Finley_Mesh* mesh=m_finleyMesh.get();
714 escriptDataC _normal=normal.getDataC();
715 switch(normal.getFunctionSpace().getTypeCode()) {
716 case(Nodes):
717 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
718 break;
719 case(ReducedNodes):
720 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
721 break;
722 case(Elements):
723 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
724 break;
725 case(ReducedElements):
726 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
727 break;
728 case (FaceElements):
729 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
730 break;
731 case (ReducedFaceElements):
732 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
733 break;
734 case(Points):
735 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
736 break;
737 case (ContactElementsOne):
738 case (ContactElementsZero):
739 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
740 break;
741 case (ReducedContactElementsOne):
742 case (ReducedContactElementsZero):
743 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
744 break;
745 case(DegreesOfFreedom):
746 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
747 break;
748 case(ReducedDegreesOfFreedom):
749 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
750 break;
751 default:
752 stringstream temp;
753 temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
754 throw FinleyAdapterException(temp.str());
755 break;
756 }
757 checkFinleyError();
758 }
759
760 //
761 // interpolates data to other domain:
762 //
763 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
764 {
765 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
766 if (targetDomain!=*this)
767 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
768
769 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
770 }
771
772 //
773 // calculates the integral of a function defined of arg:
774 //
775 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
776 {
777 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
778 if (argDomain!=*this)
779 throw FinleyAdapterException("Error - Illegal domain of integration kernel");
780
781 double blocktimer_start = blocktimer_time();
782 Finley_Mesh* mesh=m_finleyMesh.get();
783 escriptDataC _temp;
784 escript::Data temp;
785 escriptDataC _arg=arg.getDataC();
786 switch(arg.getFunctionSpace().getTypeCode()) {
787 case(Nodes):
788 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
789 _temp=temp.getDataC();
790 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
791 break;
792 case(ReducedNodes):
793 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
794 _temp=temp.getDataC();
795 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
796 break;
797 case(Elements):
798 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
799 break;
800 case(ReducedElements):
801 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
802 break;
803 case(FaceElements):
804 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
805 break;
806 case(ReducedFaceElements):
807 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
808 break;
809 case(Points):
810 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
811 break;
812 case(ContactElementsZero):
813 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
814 break;
815 case(ReducedContactElementsZero):
816 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
817 break;
818 case(ContactElementsOne):
819 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
820 break;
821 case(ReducedContactElementsOne):
822 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
823 break;
824 case(DegreesOfFreedom):
825 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
826 _temp=temp.getDataC();
827 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
828 break;
829 case(ReducedDegreesOfFreedom):
830 temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
831 _temp=temp.getDataC();
832 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
833 break;
834 default:
835 stringstream temp;
836 temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
837 throw FinleyAdapterException(temp.str());
838 break;
839 }
840 checkFinleyError();
841 blocktimer_increment("integrate()", blocktimer_start);
842 }
843
844 //
845 // calculates the gradient of arg:
846 //
847 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
848 {
849 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
850 if (argDomain!=*this)
851 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
852 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
853 if (gradDomain!=*this)
854 throw FinleyAdapterException("Error - Illegal domain of gradient");
855
856 Finley_Mesh* mesh=m_finleyMesh.get();
857 escriptDataC _grad=grad.getDataC();
858 escriptDataC nodeDataC;
859 escript::Data temp;
860 if (getMPISize()>1) {
861 if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
862 temp=escript::Data( arg, continuousFunction(asAbstractContinuousDomain()) );
863 nodeDataC = temp.getDataC();
864 } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
865 temp=escript::Data( arg, reducedContinuousFunction(asAbstractContinuousDomain()) );
866 nodeDataC = temp.getDataC();
867 } else {
868 nodeDataC = arg.getDataC();
869 }
870 } else {
871 nodeDataC = arg.getDataC();
872 }
873 switch(grad.getFunctionSpace().getTypeCode()) {
874 case(Nodes):
875 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
876 break;
877 case(ReducedNodes):
878 throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
879 break;
880 case(Elements):
881 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
882 break;
883 case(ReducedElements):
884 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
885 break;
886 case(FaceElements):
887 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
888 break;
889 case(ReducedFaceElements):
890 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
891 break;
892 case(Points):
893 throw FinleyAdapterException("Error - Gradient at points is not supported.");
894 break;
895 case(ContactElementsZero):
896 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
897 break;
898 case(ReducedContactElementsZero):
899 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
900 break;
901 case(ContactElementsOne):
902 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
903 break;
904 case(ReducedContactElementsOne):
905 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
906 break;
907 case(DegreesOfFreedom):
908 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
909 break;
910 case(ReducedDegreesOfFreedom):
911 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
912 break;
913 default:
914 stringstream temp;
915 temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
916 throw FinleyAdapterException(temp.str());
917 break;
918 }
919 checkFinleyError();
920 }
921
922 //
923 // returns the size of elements:
924 //
925 void MeshAdapter::setToSize(escript::Data& size) const
926 {
927 Finley_Mesh* mesh=m_finleyMesh.get();
928 escriptDataC tmp=size.getDataC();
929 switch(size.getFunctionSpace().getTypeCode()) {
930 case(Nodes):
931 throw FinleyAdapterException("Error - Size of nodes is not supported.");
932 break;
933 case(ReducedNodes):
934 throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
935 break;
936 case(Elements):
937 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
938 break;
939 case(ReducedElements):
940 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
941 break;
942 case(FaceElements):
943 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
944 break;
945 case(ReducedFaceElements):
946 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
947 break;
948 case(Points):
949 throw FinleyAdapterException("Error - Size of point elements is not supported.");
950 break;
951 case(ContactElementsZero):
952 case(ContactElementsOne):
953 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
954 break;
955 case(ReducedContactElementsZero):
956 case(ReducedContactElementsOne):
957 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
958 break;
959 case(DegreesOfFreedom):
960 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
961 break;
962 case(ReducedDegreesOfFreedom):
963 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
964 break;
965 default:
966 stringstream temp;
967 temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
968 throw FinleyAdapterException(temp.str());
969 break;
970 }
971 checkFinleyError();
972 }
973
974 // sets the location of nodes:
975 void MeshAdapter::setNewX(const escript::Data& new_x)
976 {
977 Finley_Mesh* mesh=m_finleyMesh.get();
978 escriptDataC tmp;
979 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
980 if (newDomain!=*this)
981 throw FinleyAdapterException("Error - Illegal domain of new point locations");
982 tmp = new_x.getDataC();
983 Finley_Mesh_setCoordinates(mesh,&tmp);
984 checkFinleyError();
985 }
986
987 // saves a data array in openDX format:
988 void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
989 {
990 int MAX_namelength=256;
991 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
992 /* win32 refactor */
993 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
994 for(int i=0;i<num_data;i++)
995 {
996 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
997 }
998
999 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1000 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
1001 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
1002
1003 boost::python::list keys=arg.keys();
1004 for (int i=0;i<num_data;++i) {
1005 std::string n=boost::python::extract<std::string>(keys[i]);
1006 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1007 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1008 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
1009 data[i]=d.getDataC();
1010 ptr_data[i]=&(data[i]);
1011 c_names[i]=&(names[i][0]);
1012 if (n.length()>MAX_namelength-1) {
1013 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
1014 c_names[i][MAX_namelength-1]='\0';
1015 } else {
1016 strcpy(c_names[i],n.c_str());
1017 }
1018 }
1019 Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1020 checkFinleyError();
1021
1022 /* win32 refactor */
1023 TMPMEMFREE(c_names);
1024 TMPMEMFREE(data);
1025 TMPMEMFREE(ptr_data);
1026 for(int i=0;i<num_data;i++)
1027 {
1028 TMPMEMFREE(names[i]);
1029 }
1030 TMPMEMFREE(names);
1031
1032 return;
1033 }
1034
1035 // saves a data array in openVTK format:
1036 void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
1037 {
1038 int MAX_namelength=256;
1039 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1040 /* win32 refactor */
1041 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1042 for(int i=0;i<num_data;i++)
1043 {
1044 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
1045 }
1046
1047 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
1048 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
1049 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
1050
1051 boost::python::list keys=arg.keys();
1052 for (int i=0;i<num_data;++i) {
1053 std::string n=boost::python::extract<std::string>(keys[i]);
1054 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1055 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
1056 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
1057 data[i]=d.getDataC();
1058 ptr_data[i]=&(data[i]);
1059 c_names[i]=&(names[i][0]);
1060 if (n.length()>MAX_namelength-1) {
1061 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
1062 c_names[i][MAX_namelength-1]='\0';
1063 } else {
1064 strcpy(c_names[i],n.c_str());
1065 }
1066 }
1067 Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1068
1069 checkFinleyError();
1070 /* win32 refactor */
1071 TMPMEMFREE(c_names);
1072 TMPMEMFREE(data);
1073 TMPMEMFREE(ptr_data);
1074 for(int i=0;i<num_data;i++)
1075 {
1076 TMPMEMFREE(names[i]);
1077 }
1078 TMPMEMFREE(names);
1079
1080 return;
1081 }
1082
1083
1084 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
1085 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1086 const int row_blocksize,
1087 const escript::FunctionSpace& row_functionspace,
1088 const int column_blocksize,
1089 const escript::FunctionSpace& column_functionspace,
1090 const int type) const
1091 {
1092 int reduceRowOrder=0;
1093 int reduceColOrder=0;
1094 // is the domain right?
1095 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
1096 if (row_domain!=*this)
1097 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1098 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
1099 if (col_domain!=*this)
1100 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1101 // is the function space type right
1102 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1103 reduceRowOrder=0;
1104 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1105 reduceRowOrder=1;
1106 } else {
1107 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1108 }
1109 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1110 reduceColOrder=0;
1111 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1112 reduceColOrder=1;
1113 } else {
1114 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1115 }
1116 // generate matrix:
1117
1118 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1119 checkFinleyError();
1120 Paso_SystemMatrix* fsystemMatrix;
1121 int trilinos = 0;
1122 if (trilinos) {
1123 #ifdef TRILINOS
1124 /* Allocation Epetra_VrbMatrix here */
1125 #endif
1126 }
1127 else {
1128 fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
1129 }
1130 checkPasoError();
1131 Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1132 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1133 }
1134
1135 //
1136 // vtkObject MeshAdapter::createVtkObject() const
1137 // TODO:
1138 //
1139
1140 //
1141 // returns true if data at the atom_type is considered as being cell centered:
1142 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1143 {
1144 switch(functionSpaceCode) {
1145 case(Nodes):
1146 case(DegreesOfFreedom):
1147 case(ReducedDegreesOfFreedom):
1148 return false;
1149 break;
1150 case(Elements):
1151 case(FaceElements):
1152 case(Points):
1153 case(ContactElementsZero):
1154 case(ContactElementsOne):
1155 case(ReducedElements):
1156 case(ReducedFaceElements):
1157 case(ReducedContactElementsZero):
1158 case(ReducedContactElementsOne):
1159 return true;
1160 break;
1161 default:
1162 stringstream temp;
1163 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1164 throw FinleyAdapterException(temp.str());
1165 break;
1166 }
1167 checkFinleyError();
1168 return false;
1169 }
1170
1171 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1172 {
1173 switch(functionSpaceType_source) {
1174 case(Nodes):
1175 switch(functionSpaceType_target) {
1176 case(Nodes):
1177 case(ReducedNodes):
1178 case(ReducedDegreesOfFreedom):
1179 case(DegreesOfFreedom):
1180 case(Elements):
1181 case(ReducedElements):
1182 case(FaceElements):
1183 case(ReducedFaceElements):
1184 case(Points):
1185 case(ContactElementsZero):
1186 case(ReducedContactElementsZero):
1187 case(ContactElementsOne):
1188 case(ReducedContactElementsOne):
1189 return true;
1190 default:
1191 stringstream temp;
1192 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1193 throw FinleyAdapterException(temp.str());
1194 }
1195 break;
1196 case(ReducedNodes):
1197 switch(functionSpaceType_target) {
1198 case(ReducedNodes):
1199 case(ReducedDegreesOfFreedom):
1200 case(Elements):
1201 case(ReducedElements):
1202 case(FaceElements):
1203 case(ReducedFaceElements):
1204 case(Points):
1205 case(ContactElementsZero):
1206 case(ReducedContactElementsZero):
1207 case(ContactElementsOne):
1208 case(ReducedContactElementsOne):
1209 return true;
1210 case(Nodes):
1211 case(DegreesOfFreedom):
1212 return false;
1213 default:
1214 stringstream temp;
1215 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1216 throw FinleyAdapterException(temp.str());
1217 }
1218 break;
1219 case(Elements):
1220 if (functionSpaceType_target==Elements) {
1221 return true;
1222 } else if (functionSpaceType_target==ReducedElements) {
1223 return true;
1224 } else {
1225 return false;
1226 }
1227 case(ReducedElements):
1228 if (functionSpaceType_target==ReducedElements) {
1229 return true;
1230 } else {
1231 return false;
1232 }
1233 case(FaceElements):
1234 if (functionSpaceType_target==FaceElements) {
1235 return true;
1236 } else if (functionSpaceType_target==ReducedFaceElements) {
1237 return true;
1238 } else {
1239 return false;
1240 }
1241 case(ReducedFaceElements):
1242 if (functionSpaceType_target==ReducedFaceElements) {
1243 return true;
1244 } else {
1245 return false;
1246 }
1247 case(Points):
1248 if (functionSpaceType_target==Points) {
1249 return true;
1250 } else {
1251 return false;
1252 }
1253 case(ContactElementsZero):
1254 case(ContactElementsOne):
1255 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1256 return true;
1257 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1258 return true;
1259 } else {
1260 return false;
1261 }
1262 case(ReducedContactElementsZero):
1263 case(ReducedContactElementsOne):
1264 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1265 return true;
1266 } else {
1267 return false;
1268 }
1269 case(DegreesOfFreedom):
1270 switch(functionSpaceType_target) {
1271 case(ReducedDegreesOfFreedom):
1272 case(DegreesOfFreedom):
1273 case(Nodes):
1274 case(ReducedNodes):
1275 case(Elements):
1276 case(ReducedElements):
1277 case(Points):
1278 case(FaceElements):
1279 case(ReducedFaceElements):
1280 case(ContactElementsZero):
1281 case(ReducedContactElementsZero):
1282 case(ContactElementsOne):
1283 case(ReducedContactElementsOne):
1284 return true;
1285 default:
1286 stringstream temp;
1287 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1288 throw FinleyAdapterException(temp.str());
1289 }
1290 break;
1291 case(ReducedDegreesOfFreedom):
1292 switch(functionSpaceType_target) {
1293 case(ReducedDegreesOfFreedom):
1294 case(ReducedNodes):
1295 case(Elements):
1296 case(ReducedElements):
1297 case(FaceElements):
1298 case(ReducedFaceElements):
1299 case(Points):
1300 case(ContactElementsZero):
1301 case(ReducedContactElementsZero):
1302 case(ContactElementsOne):
1303 case(ReducedContactElementsOne):
1304 return true;
1305 case(Nodes):
1306 case(DegreesOfFreedom):
1307 return false;
1308 default:
1309 stringstream temp;
1310 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1311 throw FinleyAdapterException(temp.str());
1312 }
1313 break;
1314 default:
1315 stringstream temp;
1316 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1317 throw FinleyAdapterException(temp.str());
1318 break;
1319 }
1320 checkFinleyError();
1321 return false;
1322 }
1323
1324 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1325 {
1326 return false;
1327 }
1328
1329 bool MeshAdapter::operator==(const AbstractDomain& other) const
1330 {
1331 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1332 if (temp!=0) {
1333 return (m_finleyMesh==temp->m_finleyMesh);
1334 } else {
1335 return false;
1336 }
1337 }
1338
1339 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1340 {
1341 return !(operator==(other));
1342 }
1343
1344 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1345 {
1346 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1347 checkPasoError();
1348 return out;
1349 }
1350
1351 escript::Data MeshAdapter::getX() const
1352 {
1353 return continuousFunction(asAbstractContinuousDomain()).getX();
1354 }
1355
1356 escript::Data MeshAdapter::getNormal() const
1357 {
1358 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1359 }
1360
1361 escript::Data MeshAdapter::getSize() const
1362 {
1363 return function(asAbstractContinuousDomain()).getSize();
1364 }
1365
1366 int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1367 {
1368 int *out=0,i;
1369 Finley_Mesh* mesh=m_finleyMesh.get();
1370 switch (functionSpaceType) {
1371 case(Nodes):
1372 out=mesh->Nodes->Id;
1373 break;
1374 case(ReducedNodes):
1375 out=mesh->Nodes->reducedNodesId;
1376 break;
1377 case(Elements):
1378 out=mesh->Elements->Id;
1379 break;
1380 case(ReducedElements):
1381 out=mesh->Elements->Id;
1382 break;
1383 case(FaceElements):
1384 out=mesh->FaceElements->Id;
1385 break;
1386 case(ReducedFaceElements):
1387 out=mesh->FaceElements->Id;
1388 break;
1389 case(Points):
1390 out=mesh->Points->Id;
1391 break;
1392 case(ContactElementsZero):
1393 out=mesh->ContactElements->Id;
1394 break;
1395 case(ReducedContactElementsZero):
1396 out=mesh->ContactElements->Id;
1397 break;
1398 case(ContactElementsOne):
1399 out=mesh->ContactElements->Id;
1400 break;
1401 case(ReducedContactElementsOne):
1402 out=mesh->ContactElements->Id;
1403 break;
1404 case(DegreesOfFreedom):
1405 out=mesh->Nodes->degreesOfFreedomId;
1406 break;
1407 case(ReducedDegreesOfFreedom):
1408 out=mesh->Nodes->reducedDegreesOfFreedomId;
1409 break;
1410 default:
1411 stringstream temp;
1412 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1413 throw FinleyAdapterException(temp.str());
1414 break;
1415 }
1416 return out;
1417 }
1418 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1419 {
1420 int out=0;
1421 Finley_Mesh* mesh=m_finleyMesh.get();
1422 switch (functionSpaceType) {
1423 case(Nodes):
1424 out=mesh->Nodes->Tag[sampleNo];
1425 break;
1426 case(ReducedNodes):
1427 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1428 break;
1429 case(Elements):
1430 out=mesh->Elements->Tag[sampleNo];
1431 break;
1432 case(ReducedElements):
1433 out=mesh->Elements->Tag[sampleNo];
1434 break;
1435 case(FaceElements):
1436 out=mesh->FaceElements->Tag[sampleNo];
1437 break;
1438 case(ReducedFaceElements):
1439 out=mesh->FaceElements->Tag[sampleNo];
1440 break;
1441 case(Points):
1442 out=mesh->Points->Tag[sampleNo];
1443 break;
1444 case(ContactElementsZero):
1445 out=mesh->ContactElements->Tag[sampleNo];
1446 break;
1447 case(ReducedContactElementsZero):
1448 out=mesh->ContactElements->Tag[sampleNo];
1449 break;
1450 case(ContactElementsOne):
1451 out=mesh->ContactElements->Tag[sampleNo];
1452 break;
1453 case(ReducedContactElementsOne):
1454 out=mesh->ContactElements->Tag[sampleNo];
1455 break;
1456 case(DegreesOfFreedom):
1457 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1458 break;
1459 case(ReducedDegreesOfFreedom):
1460 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1461 break;
1462 default:
1463 stringstream temp;
1464 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1465 throw FinleyAdapterException(temp.str());
1466 break;
1467 }
1468 return out;
1469 }
1470
1471
1472 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1473 {
1474 Finley_Mesh* mesh=m_finleyMesh.get();
1475 escriptDataC tmp=mask.getDataC();
1476 switch(functionSpaceType) {
1477 case(Nodes):
1478 Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1479 break;
1480 case(ReducedNodes):
1481 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
1482 break;
1483 case(DegreesOfFreedom):
1484 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1485 break;
1486 case(ReducedDegreesOfFreedom):
1487 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1488 break;
1489 case(Elements):
1490 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1491 break;
1492 case(ReducedElements):
1493 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1494 break;
1495 case(FaceElements):
1496 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1497 break;
1498 case(ReducedFaceElements):
1499 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1500 break;
1501 case(Points):
1502 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1503 break;
1504 case(ContactElementsZero):
1505 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1506 break;
1507 case(ReducedContactElementsZero):
1508 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1509 break;
1510 case(ContactElementsOne):
1511 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1512 break;
1513 case(ReducedContactElementsOne):
1514 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1515 break;
1516 default:
1517 stringstream temp;
1518 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1519 throw FinleyAdapterException(temp.str());
1520 }
1521 checkFinleyError();
1522 return;
1523 }
1524
1525 void MeshAdapter::setTagMap(const std::string& name, int tag)
1526 {
1527 Finley_Mesh* mesh=m_finleyMesh.get();
1528 Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
1529 checkPasoError();
1530 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1531 }
1532
1533 int MeshAdapter::getTag(const std::string& name) const
1534 {
1535 Finley_Mesh* mesh=m_finleyMesh.get();
1536 int tag=0;
1537 tag=Finley_Mesh_getTag(mesh, name.c_str());
1538 checkPasoError();
1539 // throwStandardException("MeshAdapter::getTag is not implemented.");
1540 return tag;
1541 }
1542
1543 bool MeshAdapter::isValidTagName(const std::string& name) const
1544 {
1545 Finley_Mesh* mesh=m_finleyMesh.get();
1546 return Finley_Mesh_isValidTagName(mesh,name.c_str());
1547 }
1548
1549 std::string MeshAdapter::showTagNames() const
1550 {
1551 stringstream temp;
1552 Finley_Mesh* mesh=m_finleyMesh.get();
1553 Finley_TagMap* tag_map=mesh->TagMap;
1554 while (tag_map) {
1555 temp << tag_map->name;
1556 tag_map=tag_map->next;
1557 if (tag_map) temp << ", ";
1558 }
1559 return temp.str();
1560 }
1561
1562 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26