/[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 1326 - (show annotations)
Mon Oct 1 08:10:41 2007 UTC (12 years ago) by ksteube
File size: 58844 byte(s)
Implemented domain.print_mesh_info() so we can see the distribution of elements & nodes.
Implemented -DBOUNDS_CHECK to catch an error with periodicN=True.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26