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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 210 - (show annotations)
Wed Nov 23 09:54:02 2005 UTC (14 years, 11 months ago) by robwdcock
File size: 41732 byte(s)
PARTIAL WIN32 BUILD SYSTEM AND PORT
+ All libraries now build under new build system. No unit test routines yet
+ Reversed some incorrect refactorings in Bruce.cpp

1 // $Id$
2 /*
3 ******************************************************************************
4 * *
5 * COPYRIGHT ACcESS 2004 - All Rights Reserved *
6 * *
7 * This software is the property of ACcESS. No part of this code *
8 * may be copied in any form or by any means without the expressed written *
9 * consent of ACcESS. Copying, use or modification of this software *
10 * by any unauthorised person is illegal unless that person has a software *
11 * license agreement with ACcESS. *
12 * *
13 ******************************************************************************
14 */
15 extern "C" {
16 #include "finley/finleyC/Finley.h"
17 #include "finley/finleyC/Assemble.h"
18 #include "finley/finleyC/Mesh.h"
19 #include "finley/finleyC/Finley.h"
20 #include "paso/SystemMatrix.h"
21 }
22 #include "finley/CPPAdapter/SystemMatrixAdapter.h"
23 #include "finley/CPPAdapter/MeshAdapter.h"
24 #include "finley/CPPAdapter/FinleyError.h"
25 #include "finley/CPPAdapter/FinleyAdapterException.h"
26 #include "escript/Data/FunctionSpaceFactory.h"
27 #include "escript/Data/Data.h"
28 #include "escript/Data/DataArrayView.h"
29 #include "escript/Data/FunctionSpace.h"
30 #include "escript/Data/DataFactory.h"
31 #include <boost/python/extract.hpp>
32
33 #include <iostream>
34 #include <sstream>
35
36 using namespace std;
37 using namespace escript;
38
39 namespace finley {
40
41 //
42 // define the static constants
43 MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
44 const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
45 const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
46 const int MeshAdapter::Nodes=FINLEY_NODES;
47 const int MeshAdapter::Elements=FINLEY_ELEMENTS;
48 const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
49 const int MeshAdapter::Points=FINLEY_POINTS;
50 const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
51 const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
52
53 MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
54 {
55 setFunctionSpaceTypeNames();
56 //
57 // need to use a null_deleter as Finley_Mesh_dealloc deletes the pointer
58 // for us.
59 m_finleyMesh.reset(finleyMesh,null_deleter());
60 }
61
62 //
63 // The copy constructor should just increment the use count
64 MeshAdapter::MeshAdapter(const MeshAdapter& in):
65 m_finleyMesh(in.m_finleyMesh)
66 {
67 setFunctionSpaceTypeNames();
68 }
69
70 MeshAdapter::~MeshAdapter()
71 {
72 //
73 // I hope the case for the pointer being zero has been taken care of.
74 // cout << "In MeshAdapter destructor." << endl;
75 if (m_finleyMesh.unique()) {
76 // cout << "Calling dealloc." << endl;
77 Finley_Mesh_dealloc(m_finleyMesh.get());
78 // cout << "Finished dealloc." << endl;
79 }
80 }
81
82 Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
83 return m_finleyMesh.get();
84 }
85
86 void MeshAdapter::write(const std::string& fileName) const
87 {
88 char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
89 strcpy(fName,fileName.c_str());
90 Finley_Mesh_write(m_finleyMesh.get(),fName);
91 checkFinleyError();
92 TMPMEMFREE(fName);
93 }
94
95 // void MeshAdapter::getTagList(int functionSpaceType,
96 // int* numTags) const
97 // {
98 // Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);
99 // return;
100 // }
101
102 string MeshAdapter::getDescription() const
103 {
104 return "FinleyMesh";
105 }
106
107 string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
108 {
109 FunctionSpaceNamesMapType::iterator loc;
110 loc=m_functionSpaceTypeNames.find(functionSpaceType);
111 if (loc==m_functionSpaceTypeNames.end()) {
112 return "Invalid function space type code.";
113 } else {
114 return loc->second;
115 }
116 }
117
118 bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
119 {
120 FunctionSpaceNamesMapType::iterator loc;
121 loc=m_functionSpaceTypeNames.find(functionSpaceType);
122 return (loc!=m_functionSpaceTypeNames.end());
123 }
124
125 void MeshAdapter::setFunctionSpaceTypeNames()
126 {
127 m_functionSpaceTypeNames.insert
128 (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));
129 m_functionSpaceTypeNames.insert
130 (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));
131 m_functionSpaceTypeNames.insert
132 (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
133 m_functionSpaceTypeNames.insert
134 (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
135 m_functionSpaceTypeNames.insert
136 (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
137 m_functionSpaceTypeNames.insert
138 (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
139 m_functionSpaceTypeNames.insert
140 (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
141 m_functionSpaceTypeNames.insert
142 (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
143 }
144
145 int MeshAdapter::getContinuousFunctionCode() const
146 {
147 return Nodes;
148 }
149
150 int MeshAdapter::getFunctionCode() const
151 {
152 return Elements;
153 }
154
155 int MeshAdapter::getFunctionOnBoundaryCode() const
156 {
157 return FaceElements;
158 }
159
160 int MeshAdapter::getFunctionOnContactZeroCode() const
161 {
162 return ContactElementsZero;
163 }
164
165 int MeshAdapter::getFunctionOnContactOneCode() const
166 {
167 return ContactElementsOne;
168 }
169
170 int MeshAdapter::getSolutionCode() const
171 {
172 return DegreesOfFreedom;
173 }
174
175 int MeshAdapter::getReducedSolutionCode() const
176 {
177 return ReducedDegreesOfFreedom;
178 }
179
180 int MeshAdapter::getDiracDeltaFunctionCode() const
181 {
182 return Points;
183 }
184
185 //
186 // returns a pointer to the tag list of samples of functionSpaceType
187 //
188 void MeshAdapter::getTagList(int functionSpaceType, int** tagList,
189 int* numTags) const
190 {
191 *tagList=NULL;
192 *numTags=0;
193 Finley_Mesh* mesh=m_finleyMesh.get();
194 switch (functionSpaceType) {
195 case(Nodes):
196 if (mesh->Nodes!=NULL) {
197 *tagList=mesh->Nodes->Tag;
198 *numTags=mesh->Nodes->numNodes;
199 }
200 break;
201 case(Elements):
202 if (mesh->Elements!=NULL) {
203 *tagList=mesh->Elements->Tag;
204 *numTags=mesh->Elements->numElements;
205 }
206 break;
207 case(FaceElements):
208 if (mesh->FaceElements!=NULL) {
209 *tagList=mesh->FaceElements->Tag;
210 *numTags=mesh->FaceElements->numElements;
211 }
212 break;
213 case(Points):
214 if (mesh->Points!=NULL) {
215 *tagList=mesh->Points->Tag;
216 *numTags=mesh->Points->numElements;
217 }
218 break;
219 case(ContactElementsZero):
220 if (mesh->ContactElements!=NULL) {
221 *tagList=mesh->ContactElements->Tag;
222 *numTags=mesh->ContactElements->numElements;
223 }
224 break;
225 case(ContactElementsOne):
226 if (mesh->ContactElements!=NULL) {
227 *tagList=mesh->ContactElements->Tag;
228 *numTags=mesh->ContactElements->numElements;
229 }
230 break;
231 case(DegreesOfFreedom):
232 if (mesh->Nodes!=NULL) {
233 *tagList=NULL;
234 *numTags=0;
235 }
236 break;
237 case(ReducedDegreesOfFreedom):
238 if (mesh->Nodes!=NULL) {
239 *tagList=NULL;
240 *numTags=0;
241 }
242 break;
243 default:
244 stringstream temp;
245 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
246 throw FinleyAdapterException(temp.str());
247 break;
248 }
249 if (*tagList==NULL) {
250 stringstream temp;
251 temp << "Error - no tags available for " << functionSpaceType << " for domain: " << getDescription();
252 throw FinleyAdapterException(temp.str());
253 }
254 return;
255 }
256
257 //
258 // returns a pointer to the reference no list of samples of functionSpaceType
259 //
260 void MeshAdapter::getReferenceNoList(int functionSpaceType, int** referenceNoList,
261 int* numReferenceNo) const
262 {
263 *referenceNoList=NULL;
264 *numReferenceNo=0;
265 Finley_Mesh* mesh=m_finleyMesh.get();
266 switch (functionSpaceType) {
267 case(Nodes):
268 if (mesh->Nodes!=NULL) {
269 *referenceNoList=mesh->Nodes->Id;
270 *numReferenceNo=mesh->Nodes->numNodes;
271 }
272 break;
273 case(Elements):
274 if (mesh->Elements!=NULL) {
275 *referenceNoList=mesh->Elements->Id;
276 *numReferenceNo=mesh->Elements->numElements;
277 }
278 break;
279 case(FaceElements):
280 if (mesh->FaceElements!=NULL) {
281 *referenceNoList=mesh->FaceElements->Id;
282 *numReferenceNo=mesh->FaceElements->numElements;
283 }
284 break;
285 case(Points):
286 if (mesh->Points!=NULL) {
287 *referenceNoList=mesh->Points->Id;
288 *numReferenceNo=mesh->Points->numElements;
289 }
290 break;
291 case(ContactElementsZero):
292 if (mesh->ContactElements!=NULL) {
293 *referenceNoList=mesh->ContactElements->Id;
294 *numReferenceNo=mesh->ContactElements->numElements;
295 }
296 break;
297 case(ContactElementsOne):
298 if (mesh->ContactElements!=NULL) {
299 *referenceNoList=mesh->ContactElements->Id;
300 *numReferenceNo=mesh->ContactElements->numElements;
301 }
302 break;
303 case(DegreesOfFreedom):
304 if (mesh->Nodes!=NULL) {
305 *referenceNoList=NULL;
306 *numReferenceNo=0;
307 }
308 break;
309 case(ReducedDegreesOfFreedom):
310 if (mesh->Nodes!=NULL) {
311 *referenceNoList=NULL;
312 *numReferenceNo=0;
313 }
314 break;
315 default:
316 stringstream temp;
317 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
318 throw FinleyAdapterException(temp.str());
319 break;
320 }
321 if (*referenceNoList==NULL) {
322 stringstream temp;
323 temp << "Error - reference number list available for " << functionSpaceType << " for domain: " << getDescription();
324 throw FinleyAdapterException(temp.str());
325 }
326 return;
327 }
328
329 //
330 // return the spatial dimension of the Mesh:
331 //
332 int MeshAdapter::getDim() const
333 {
334 int numDim=Finley_Mesh_getDim(m_finleyMesh.get());
335 checkFinleyError();
336 return numDim;
337 }
338
339 //
340 // return the number of data points per sample and the number of samples
341 // needed to represent data on a parts of the mesh.
342 //
343 pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
344 {
345 int numDataPointsPerSample=0;
346 int numSamples=0;
347 Finley_Mesh* mesh=m_finleyMesh.get();
348 switch (functionSpaceCode) {
349 case(Nodes):
350 numDataPointsPerSample=1;
351 if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
352 break;
353 case(Elements):
354 if (mesh->Elements!=NULL) {
355 numSamples=mesh->Elements->numElements;
356 numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
357 }
358 break;
359 case(FaceElements):
360 if (mesh->FaceElements!=NULL) {
361 numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
362 numSamples=mesh->FaceElements->numElements;
363 }
364 break;
365 case(Points):
366 if (mesh->Points!=NULL) {
367 numDataPointsPerSample=1;
368 numSamples=mesh->Points->numElements;
369 }
370 break;
371 case(ContactElementsZero):
372 if (mesh->ContactElements!=NULL) {
373 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
374 numSamples=mesh->ContactElements->numElements;
375 }
376 break;
377 case(ContactElementsOne):
378 if (mesh->ContactElements!=NULL) {
379 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
380 numSamples=mesh->ContactElements->numElements;
381 }
382 break;
383 case(DegreesOfFreedom):
384 if (mesh->Nodes!=NULL) {
385 numDataPointsPerSample=1;
386 numSamples=mesh->Nodes->numDegreesOfFreedom;
387 }
388 break;
389 case(ReducedDegreesOfFreedom):
390 if (mesh->Nodes!=NULL) {
391 numDataPointsPerSample=1;
392 numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
393 }
394 break;
395 default:
396 stringstream temp;
397 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
398 throw FinleyAdapterException(temp.str());
399 break;
400 }
401 return pair<int,int>(numDataPointsPerSample,numSamples);
402 }
403
404 //
405 // adds linear PDE of second order into a given stiffness matrix and right hand side:
406 //
407 void MeshAdapter::addPDEToSystem(
408 SystemMatrixAdapter& mat, Data& rhs,
409 const Data& A, const Data& B, const Data& C,const Data& D,const Data& X,const Data& Y,
410 const Data& d, const Data& y,
411 const Data& d_contact,const Data& y_contact) const
412 {
413 Finley_Mesh* mesh=m_finleyMesh.get();
414 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(),&(rhs.getDataC()),
415 &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));
416 checkFinleyError();
417 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,
418 mat.getPaso_SystemMatrix(),
419 &(rhs.getDataC()),
420 &(d.getDataC()),&(y.getDataC()),
421 Finley_Assemble_handelShapeMissMatch_Mean_out);
422 checkFinleyError();
423 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,
424 mat.getPaso_SystemMatrix(),
425 &(rhs.getDataC()),
426 &(d_contact.getDataC()),
427 &(y_contact.getDataC()),
428 Finley_Assemble_handelShapeMissMatch_Step_out);
429 checkFinleyError();
430 }
431
432 //
433 // adds linear PDE of second order into the right hand side only
434 //
435 void MeshAdapter::addPDEToRHS( Data& rhs,
436 const Data& X,const Data& Y, const Data& y, const Data& y_contact) const
437 {
438 Finley_Mesh* mesh=m_finleyMesh.get();
439
440 // Finley_Assemble_PDE_RHS(mesh->Nodes,mesh->Elements,&(rhs.getDataC()),&(X.getDataC()),&(Y.getDataC()));
441 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));
442 checkFinleyError();
443
444 // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->FaceElements,&(rhs.getDataC()),&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);
445 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);
446
447 checkFinleyError();
448 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);
449 // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->ContactElements,&(rhs.getDataC()),&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);
450 checkFinleyError();
451 }
452
453 //
454 // interpolates data between different function spaces:
455 //
456 void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const
457 {
458 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
459 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
460 if (inDomain!=*this)
461 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
462 if (targetDomain!=*this)
463 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
464
465 Finley_Mesh* mesh=m_finleyMesh.get();
466 switch(in.getFunctionSpace().getTypeCode()) {
467 case(Nodes):
468 switch(target.getFunctionSpace().getTypeCode()) {
469 case(Nodes):
470 case(ReducedDegreesOfFreedom):
471 case(DegreesOfFreedom):
472 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
473 break;
474 case(Elements):
475 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
476 break;
477 case(FaceElements):
478 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
479 break;
480 case(Points):
481 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
482 break;
483 case(ContactElementsZero):
484 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
485 break;
486 case(ContactElementsOne):
487 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
488 break;
489 default:
490 stringstream temp;
491 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
492 throw FinleyAdapterException(temp.str());
493 break;
494 }
495 break;
496 case(Elements):
497 if (target.getFunctionSpace().getTypeCode()==Elements) {
498 Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));
499 } else {
500 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
501 }
502 break;
503 case(FaceElements):
504 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
505 Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
506 } else {
507 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
508 break;
509 }
510 case(Points):
511 if (target.getFunctionSpace().getTypeCode()==Points) {
512 Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
513 } else {
514 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
515 break;
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.getDataC()),&(in.getDataC()));
522 } else {
523 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
524 break;
525 }
526 break;
527 case(DegreesOfFreedom):
528 switch(target.getFunctionSpace().getTypeCode()) {
529 case(ReducedDegreesOfFreedom):
530 case(DegreesOfFreedom):
531 case(Nodes):
532 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
533 break;
534 case(Elements):
535 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
536 break;
537 case(FaceElements):
538 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
539 break;
540 case(Points):
541 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
542 break;
543 case(ContactElementsZero):
544 case(ContactElementsOne):
545 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
546 break;
547 default:
548 stringstream temp;
549 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
550 throw FinleyAdapterException(temp.str());
551 break;
552 }
553 break;
554 case(ReducedDegreesOfFreedom):
555 switch(target.getFunctionSpace().getTypeCode()) {
556 case(ReducedDegreesOfFreedom):
557 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
558 break;
559 case(Elements):
560 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
561 break;
562 case(FaceElements):
563 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
564 break;
565 case(Points):
566 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
567 break;
568 case(ContactElementsZero):
569 case(ContactElementsOne):
570 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
571 break;
572 case(Nodes):
573 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
574 break;
575 case(DegreesOfFreedom):
576 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
577 break;
578 default:
579 stringstream temp;
580 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
581 throw FinleyAdapterException(temp.str());
582 break;
583 }
584 break;
585 default:
586 stringstream temp;
587 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
588 throw FinleyAdapterException(temp.str());
589 break;
590 }
591 checkFinleyError();
592 }
593
594 //
595 // copies the locations of sample points into x:
596 //
597 void MeshAdapter::setToX(Data& arg) const
598 {
599 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
600 if (argDomain!=*this)
601 throw FinleyAdapterException("Error - Illegal domain of data point locations");
602 Finley_Mesh* mesh=m_finleyMesh.get();
603 // in case of values node coordinates we can do the job directly:
604 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
605 Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
606 } else {
607 Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
608 Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
609 // this is then interpolated onto arg:
610 interpolateOnDomain(arg,tmp_data);
611 }
612 checkFinleyError();
613 }
614
615 //
616 // return the normal vectors at the location of data points as a Data object:
617 //
618 void MeshAdapter::setToNormal(Data& normal) const
619 {
620 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
621 if (normalDomain!=*this)
622 throw FinleyAdapterException("Error - Illegal domain of normal locations");
623 Finley_Mesh* mesh=m_finleyMesh.get();
624 switch(normal.getFunctionSpace().getTypeCode()) {
625 case(Nodes):
626 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
627 break;
628 case(Elements):
629 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
630 break;
631 case (FaceElements):
632 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
633 break;
634 case(Points):
635 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
636 break;
637 case (ContactElementsOne):
638 case (ContactElementsZero):
639 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
640 break;
641 case(DegreesOfFreedom):
642 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
643 break;
644 case(ReducedDegreesOfFreedom):
645 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
646 break;
647 default:
648 stringstream temp;
649 temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
650 throw FinleyAdapterException(temp.str());
651 break;
652 }
653 checkFinleyError();
654 }
655
656 //
657 // interpolates data to other domain:
658 //
659 void MeshAdapter::interpolateACross(Data& target,const Data& source) const
660 {
661 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
662 if (targetDomain!=*this)
663 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
664
665 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
666 }
667
668 //
669 // calculates the integral of a function defined of arg:
670 //
671 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const
672 {
673 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
674 if (argDomain!=*this)
675 throw FinleyAdapterException("Error - Illegal domain of integration kernel");
676
677 Finley_Mesh* mesh=m_finleyMesh.get();
678 switch(arg.getFunctionSpace().getTypeCode()) {
679 case(Nodes):
680 throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");
681 break;
682 case(Elements):
683 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);
684 break;
685 case(FaceElements):
686 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);
687 break;
688 case(Points):
689 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
690 break;
691 case(ContactElementsZero):
692 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
693 break;
694 case(ContactElementsOne):
695 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
696 break;
697 case(DegreesOfFreedom):
698 throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");
699 break;
700 case(ReducedDegreesOfFreedom):
701 throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");
702 break;
703 default:
704 stringstream temp;
705 temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
706 throw FinleyAdapterException(temp.str());
707 break;
708 }
709 checkFinleyError();
710 }
711
712 //
713 // calculates the gradient of arg:
714 //
715 void MeshAdapter::setToGradient(Data& grad,const Data& arg) const
716 {
717 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
718 if (argDomain!=*this)
719 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
720 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
721 if (gradDomain!=*this)
722 throw FinleyAdapterException("Error - Illegal domain of gradient");
723
724 Finley_Mesh* mesh=m_finleyMesh.get();
725 switch(grad.getFunctionSpace().getTypeCode()) {
726 case(Nodes):
727 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
728 break;
729 case(Elements):
730 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));
731 break;
732 case(FaceElements):
733 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));
734 break;
735 case(Points):
736 throw FinleyAdapterException("Error - Gradient at points is not supported.");
737 break;
738 case(ContactElementsZero):
739 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
740 break;
741 case(ContactElementsOne):
742 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
743 break;
744 case(DegreesOfFreedom):
745 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
746 break;
747 case(ReducedDegreesOfFreedom):
748 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
749 break;
750 default:
751 stringstream temp;
752 temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
753 throw FinleyAdapterException(temp.str());
754 break;
755 }
756 checkFinleyError();
757 }
758
759 //
760 // returns the size of elements:
761 //
762 void MeshAdapter::setToSize(Data& size) const
763 {
764 Finley_Mesh* mesh=m_finleyMesh.get();
765 escriptDataC tmp=size.getDataC();
766 switch(size.getFunctionSpace().getTypeCode()) {
767 case(Nodes):
768 throw FinleyAdapterException("Error - Size of nodes is not supported.");
769 break;
770 case(Elements):
771 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
772 break;
773 case(FaceElements):
774 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
775 break;
776 case(Points):
777 throw FinleyAdapterException("Error - Size of point elements is not supported.");
778 break;
779 case(ContactElementsZero):
780 case(ContactElementsOne):
781 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
782 break;
783 case(DegreesOfFreedom):
784 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
785 break;
786 case(ReducedDegreesOfFreedom):
787 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
788 break;
789 default:
790 stringstream temp;
791 temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
792 throw FinleyAdapterException(temp.str());
793 break;
794 }
795 checkFinleyError();
796 }
797
798 // sets the location of nodes:
799 void MeshAdapter::setNewX(const Data& new_x)
800 {
801 Finley_Mesh* mesh=m_finleyMesh.get();
802 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
803 if (newDomain!=*this)
804 throw FinleyAdapterException("Error - Illegal domain of new point locations");
805 Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));
806 checkFinleyError();
807 }
808
809 // saves a data array in openDX format:
810 void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
811 {
812 int MAX_namelength=256;
813 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
814 /* win32 refactor */
815 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
816 for(int i=0;i<num_data;i++)
817 {
818 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
819 }
820
821 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
822 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
823 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
824
825 boost::python::list keys=arg.keys();
826 for (int i=0;i<num_data;++i) {
827 Data& d=boost::python::extract<Data&>(arg[keys[i]]);
828 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
829 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
830 data[i]=d.getDataC();
831 ptr_data[i]=&(data[i]);
832 std::string n=boost::python::extract<std::string>(keys[i]);
833 c_names[i]=&(names[i][0]);
834 if (n.length()>MAX_namelength-1) {
835 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
836 c_names[i][MAX_namelength-1]='\0';
837 } else {
838 strcpy(c_names[i],n.c_str());
839 }
840 }
841 Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
842 checkFinleyError();
843
844 /* win32 refactor */
845 TMPMEMFREE(c_names);
846 TMPMEMFREE(data);
847 TMPMEMFREE(ptr_data);
848 for(int i=0;i<num_data;i++)
849 {
850 TMPMEMFREE(names[i]);
851 }
852 TMPMEMFREE(names);
853
854 return;
855 }
856
857 // saves a data array in openVTK format:
858 void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
859 {
860 int MAX_namelength=256;
861 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
862 /* win32 refactor */
863 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
864 for(int i=0;i<num_data;i++)
865 {
866 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
867 }
868
869 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
870 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
871 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
872
873 boost::python::list keys=arg.keys();
874 for (int i=0;i<num_data;++i) {
875 Data& d=boost::python::extract<Data&>(arg[keys[i]]);
876 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
877 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
878 data[i]=d.getDataC();
879 ptr_data[i]=&(data[i]);
880 std::string n=boost::python::extract<std::string>(keys[i]);
881 c_names[i]=&(names[i][0]);
882 if (n.length()>MAX_namelength-1) {
883 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
884 c_names[i][MAX_namelength-1]='\0';
885 } else {
886 strcpy(c_names[i],n.c_str());
887 }
888 }
889 Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
890 checkFinleyError();
891 /* win32 refactor */
892 TMPMEMFREE(c_names);
893 TMPMEMFREE(data);
894 TMPMEMFREE(ptr_data);
895 for(int i=0;i<num_data;i++)
896 {
897 TMPMEMFREE(names[i]);
898 }
899 TMPMEMFREE(names);
900
901 return;
902 }
903
904
905 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
906 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
907 const int row_blocksize,
908 const escript::FunctionSpace& row_functionspace,
909 const int column_blocksize,
910 const escript::FunctionSpace& column_functionspace,
911 const int type) const
912 {
913 int reduceRowOrder=0;
914 int reduceColOrder=0;
915 // is the domain right?
916 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
917 if (row_domain!=*this)
918 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
919 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
920 if (col_domain!=*this)
921 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
922 // is the function space type right
923 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
924 reduceRowOrder=0;
925 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
926 reduceRowOrder=1;
927 } else {
928 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
929 }
930 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
931 reduceColOrder=0;
932 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
933 reduceColOrder=1;
934 } else {
935 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
936 }
937 // generate matrix:
938
939 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
940 checkFinleyError();
941 Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
942 checkPasoError();
943 Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
944 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
945 }
946
947 //
948 // vtkObject MeshAdapter::createVtkObject() const
949 // TODO:
950 //
951
952 //
953 // returns true if data at the atom_type is considered as being cell centered:
954 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
955 {
956 switch(functionSpaceCode) {
957 case(Nodes):
958 case(DegreesOfFreedom):
959 case(ReducedDegreesOfFreedom):
960 return false;
961 break;
962 case(Elements):
963 case(FaceElements):
964 case(Points):
965 case(ContactElementsZero):
966 case(ContactElementsOne):
967 return true;
968 break;
969 default:
970 stringstream temp;
971 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
972 throw FinleyAdapterException(temp.str());
973 break;
974 }
975 checkFinleyError();
976 return false;
977 }
978
979 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
980 {
981 switch(functionSpaceType_source) {
982 case(Nodes):
983 switch(functionSpaceType_target) {
984 case(Nodes):
985 case(ReducedDegreesOfFreedom):
986 case(DegreesOfFreedom):
987 case(Elements):
988 case(FaceElements):
989 case(Points):
990 case(ContactElementsZero):
991 case(ContactElementsOne):
992 return true;
993 default:
994 stringstream temp;
995 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
996 throw FinleyAdapterException(temp.str());
997 }
998 break;
999 case(Elements):
1000 if (functionSpaceType_target==Elements) {
1001 return true;
1002 } else {
1003 return false;
1004 }
1005 case(FaceElements):
1006 if (functionSpaceType_target==FaceElements) {
1007 return true;
1008 } else {
1009 return false;
1010 }
1011 case(Points):
1012 if (functionSpaceType_target==Points) {
1013 return true;
1014 } else {
1015 return false;
1016 }
1017 case(ContactElementsZero):
1018 case(ContactElementsOne):
1019 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1020 return true;
1021 } else {
1022 return false;
1023 }
1024 case(DegreesOfFreedom):
1025 switch(functionSpaceType_target) {
1026 case(ReducedDegreesOfFreedom):
1027 case(DegreesOfFreedom):
1028 case(Nodes):
1029 case(Elements):
1030 case(FaceElements):
1031 case(Points):
1032 case(ContactElementsZero):
1033 case(ContactElementsOne):
1034 return true;
1035 default:
1036 stringstream temp;
1037 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1038 throw FinleyAdapterException(temp.str());
1039 }
1040 break;
1041 case(ReducedDegreesOfFreedom):
1042 switch(functionSpaceType_target) {
1043 case(ReducedDegreesOfFreedom):
1044 case(Elements):
1045 case(FaceElements):
1046 case(Points):
1047 case(ContactElementsZero):
1048 case(ContactElementsOne):
1049 return true;
1050 case(Nodes):
1051 case(DegreesOfFreedom):
1052 return false;
1053 default:
1054 stringstream temp;
1055 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1056 throw FinleyAdapterException(temp.str());
1057 }
1058 break;
1059 default:
1060 stringstream temp;
1061 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1062 throw FinleyAdapterException(temp.str());
1063 break;
1064 }
1065 checkFinleyError();
1066 return false;
1067 }
1068
1069 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1070 {
1071 return false;
1072 }
1073
1074 bool MeshAdapter::operator==(const AbstractDomain& other) const
1075 {
1076 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1077 if (temp!=0) {
1078 return (m_finleyMesh==temp->m_finleyMesh);
1079 } else {
1080 return false;
1081 }
1082 }
1083
1084 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1085 {
1086 return !(operator==(other));
1087 }
1088
1089 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1090 {
1091 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1092 checkPasoError();
1093 return out;
1094 }
1095
1096 Data MeshAdapter::getX() const
1097 {
1098 return continuousFunction(asAbstractContinuousDomain()).getX();
1099 }
1100
1101 Data MeshAdapter::getNormal() const
1102 {
1103 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1104 }
1105
1106 Data MeshAdapter::getSize() const
1107 {
1108 return function(asAbstractContinuousDomain()).getSize();
1109 }
1110
1111 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1112 {
1113 int* tagList;
1114 int numTags;
1115 getTagList(functionSpaceType, &tagList, &numTags);
1116 return tagList[sampleNo];
1117 }
1118
1119 int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const
1120 {
1121 int* referenceNoList;
1122 int numReferenceNo;
1123 getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);
1124 return referenceNoList[sampleNo];
1125 }
1126
1127 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26