/[escript]/trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Contents of /trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 148 - (show annotations)
Tue Aug 23 01:24:31 2005 UTC (17 years, 7 months ago) by jgs
File size: 38264 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-08-23

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26