/[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 150 - (show annotations)
Thu Sep 15 03:44:45 2005 UTC (13 years, 10 months ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 38207 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26