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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 149 - (show annotations)
Thu Sep 1 03:31:39 2005 UTC (13 years, 7 months ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 38304 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

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
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: "
245 << 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: "
318 << functionSpaceType << " for domain: " << getDescription();
319 throw FinleyAdapterException(temp.str());
320 break;
321 }
322 if (*referenceNoList==NULL) {
323 stringstream temp;
324 temp << "Error - reference number list available for " << functionSpaceType << " for domain: " << getDescription();
325 throw FinleyAdapterException(temp.str());
326 }
327 return;
328 }
329
330 //
331 // return the spatial dimension of the Mesh:
332 //
333 int MeshAdapter::getDim() const
334 {
335 int numDim=Finley_Mesh_getDim(m_finleyMesh.get());
336 checkFinleyError();
337 return numDim;
338 }
339
340 //
341 // return the number of data points per sample and the number of samples
342 // needed to represent data on a parts of the mesh.
343 //
344 pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
345 {
346 int numDataPointsPerSample=0;
347 int numSamples=0;
348 Finley_Mesh* mesh=m_finleyMesh.get();
349 switch (functionSpaceCode) {
350 case(Nodes):
351 numDataPointsPerSample=1;
352 if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
353 break;
354 case(Elements):
355 if (mesh->Elements!=NULL) {
356 numSamples=mesh->Elements->numElements;
357 numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
358 }
359 break;
360 case(FaceElements):
361 if (mesh->FaceElements!=NULL) {
362 numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
363 numSamples=mesh->FaceElements->numElements;
364 }
365 break;
366 case(Points):
367 if (mesh->Points!=NULL) {
368 numDataPointsPerSample=1;
369 numSamples=mesh->Points->numElements;
370 }
371 break;
372 case(ContactElementsZero):
373 if (mesh->ContactElements!=NULL) {
374 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
375 numSamples=mesh->ContactElements->numElements;
376 }
377 break;
378 case(ContactElementsOne):
379 if (mesh->ContactElements!=NULL) {
380 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
381 numSamples=mesh->ContactElements->numElements;
382 }
383 break;
384 case(DegreesOfFreedom):
385 if (mesh->Nodes!=NULL) {
386 numDataPointsPerSample=1;
387 numSamples=mesh->Nodes->numDegreesOfFreedom;
388 }
389 break;
390 case(ReducedDegreesOfFreedom):
391 if (mesh->Nodes!=NULL) {
392 numDataPointsPerSample=1;
393 numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
394 }
395 break;
396 default:
397 stringstream temp;
398 temp << "Error - Invalid function space type: "
399 << functionSpaceCode << " for domain: " << getDescription();
400 throw FinleyAdapterException(temp.str());
401 break;
402 }
403 return pair<int,int>(numDataPointsPerSample,numSamples);
404 }
405
406 //
407 // adds linear PDE of second order into a given stiffness matrix and right hand side:
408 //
409 void MeshAdapter::addPDEToSystem(
410 SystemMatrixAdapter& mat, Data& rhs,
411 const Data& A, const Data& B, const Data& C,const Data& D,const Data& X,const Data& Y,
412 const Data& d, const Data& y,
413 const Data& d_contact,const Data& y_contact) const
414 {
415 Finley_Mesh* mesh=m_finleyMesh.get();
416 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getFinley_SystemMatrix(),&(rhs.getDataC()),
417 &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));
418 checkFinleyError();
419 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,
420 mat.getFinley_SystemMatrix(),
421 &(rhs.getDataC()),
422 &(d.getDataC()),&(y.getDataC()),
423 Finley_Assemble_handelShapeMissMatch_Mean_out);
424 checkFinleyError();
425 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,
426 mat.getFinley_SystemMatrix(),
427 &(rhs.getDataC()),
428 &(d_contact.getDataC()),
429 &(y_contact.getDataC()),
430 Finley_Assemble_handelShapeMissMatch_Step_out);
431 checkFinleyError();
432 }
433
434 //
435 // adds linear PDE of second order into the right hand side only
436 //
437 void MeshAdapter::addPDEToRHS( Data& rhs,
438 const Data& X,const Data& Y, const Data& y, const Data& y_contact) const
439 {
440 Finley_Mesh* mesh=m_finleyMesh.get();
441
442 // Finley_Assemble_PDE_RHS(mesh->Nodes,mesh->Elements,&(rhs.getDataC()),&(X.getDataC()),&(Y.getDataC()));
443 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));
444 checkFinleyError();
445
446 // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->FaceElements,&(rhs.getDataC()),&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);
447 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);
448
449 checkFinleyError();
450 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);
451 // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->ContactElements,&(rhs.getDataC()),&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);
452 checkFinleyError();
453 }
454
455 //
456 // interpolates data between different function spaces:
457 //
458 void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const
459 {
460 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
461 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
462 if (inDomain!=*this)
463 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
464 if (targetDomain!=*this)
465 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
466
467 Finley_Mesh* mesh=m_finleyMesh.get();
468 switch(in.getFunctionSpace().getTypeCode()) {
469 case(Nodes):
470 switch(target.getFunctionSpace().getTypeCode()) {
471 case(Nodes):
472 case(ReducedDegreesOfFreedom):
473 case(DegreesOfFreedom):
474 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
475 break;
476 case(Elements):
477 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
478 break;
479 case(FaceElements):
480 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
481 break;
482 case(Points):
483 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
484 break;
485 case(ContactElementsZero):
486 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
487 break;
488 case(ContactElementsOne):
489 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
490 break;
491 default:
492 Finley_ErrorCode=TYPE_ERROR;
493 sprintf(Finley_ErrorMsg,"Interpolation on Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());
494 break;
495 }
496 break;
497 case(Elements):
498 if (target.getFunctionSpace().getTypeCode()==Elements) {
499 Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));
500 } else {
501 Finley_ErrorCode=TYPE_ERROR;
502 sprintf(Finley_ErrorMsg,"No interpolation with data on elements possible.");
503 }
504 break;
505 case(FaceElements):
506 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
507 Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
508 } else {
509 Finley_ErrorCode=TYPE_ERROR;
510 sprintf(Finley_ErrorMsg,"No interpolation with data on face elements possible.");
511 break;
512 }
513 case(Points):
514 if (target.getFunctionSpace().getTypeCode()==Points) {
515 Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
516 } else {
517 Finley_ErrorCode=TYPE_ERROR;
518 sprintf(Finley_ErrorMsg,"No interpolation with data on points possible.");
519 break;
520 }
521 break;
522 case(ContactElementsZero):
523 case(ContactElementsOne):
524 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
525 Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));
526 } else {
527 Finley_ErrorCode=TYPE_ERROR;
528 sprintf(Finley_ErrorMsg,"No interpolation with data on contact elements possible.");
529 break;
530 }
531 break;
532 case(DegreesOfFreedom):
533 switch(target.getFunctionSpace().getTypeCode()) {
534 case(ReducedDegreesOfFreedom):
535 case(DegreesOfFreedom):
536 case(Nodes):
537 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
538 break;
539 case(Elements):
540 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
541 break;
542 case(FaceElements):
543 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
544 break;
545 case(Points):
546 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
547 break;
548 case(ContactElementsZero):
549 case(ContactElementsOne):
550 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
551 break;
552 default:
553 Finley_ErrorCode=TYPE_ERROR;
554 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());
555 break;
556 }
557 break;
558 case(ReducedDegreesOfFreedom):
559 switch(target.getFunctionSpace().getTypeCode()) {
560 case(ReducedDegreesOfFreedom):
561 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
562 break;
563 case(Elements):
564 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
565 break;
566 case(FaceElements):
567 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
568 break;
569 case(Points):
570 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
571 break;
572 case(ContactElementsZero):
573 case(ContactElementsOne):
574 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
575 break;
576 case(Nodes):
577 Finley_ErrorCode=TYPE_ERROR;
578 sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
579 break;
580 case(DegreesOfFreedom):
581 Finley_ErrorCode=TYPE_ERROR;
582 sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
583 break;
584 default:
585 Finley_ErrorCode=TYPE_ERROR;
586 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());
587 break;
588 }
589 break;
590 default:
591 Finley_ErrorCode=TYPE_ERROR;
592 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",in.getFunctionSpace().getTypeCode());
593 break;
594 }
595 checkFinleyError();
596 }
597
598 //
599 // copies the locations of sample points into x:
600 //
601 void MeshAdapter::setToX(Data& arg) const
602 {
603 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
604 if (argDomain!=*this)
605 throw FinleyAdapterException("Error - Illegal domain of data point locations");
606 Finley_Mesh* mesh=m_finleyMesh.get();
607 // in case of values node coordinates we can do the job directly:
608 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
609 Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
610 } else {
611 Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
612 Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
613 // this is then interpolated onto arg:
614 interpolateOnDomain(arg,tmp_data);
615 }
616 checkFinleyError();
617 }
618
619 //
620 // return the normal vectors at the location of data points as a Data object:
621 //
622 void MeshAdapter::setToNormal(Data& normal) const
623 {
624 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
625 if (normalDomain!=*this)
626 throw FinleyAdapterException("Error - Illegal domain of normal locations");
627 Finley_Mesh* mesh=m_finleyMesh.get();
628 switch(normal.getFunctionSpace().getTypeCode()) {
629 case(Nodes):
630 Finley_ErrorCode=TYPE_ERROR;
631 sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for nodes");
632 break;
633 case(Elements):
634 Finley_ErrorCode=TYPE_ERROR;
635 sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for elements");
636 break;
637 case (FaceElements):
638 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
639 break;
640 case(Points):
641 Finley_ErrorCode=TYPE_ERROR;
642 sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for point elements");
643 break;
644 case (ContactElementsOne):
645 case (ContactElementsZero):
646 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
647 break;
648 case(DegreesOfFreedom):
649 Finley_ErrorCode=TYPE_ERROR;
650 sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for degrees of freedom.");
651 break;
652 case(ReducedDegreesOfFreedom):
653 Finley_ErrorCode=TYPE_ERROR;
654 sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for reduced degrees of freedom.");
655 break;
656 default:
657 Finley_ErrorCode=TYPE_ERROR;
658 sprintf(Finley_ErrorMsg,"Normal Vectors: Finley does not know anything about function space type %d",normal.getFunctionSpace().getTypeCode());
659 break;
660 }
661 checkFinleyError();
662 }
663
664 //
665 // interpolates data to other domain:
666 //
667 void MeshAdapter::interpolateACross(Data& target,const Data& source) const
668 {
669 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
670 if (targetDomain!=*this)
671 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
672
673 Finley_ErrorCode=SYSTEM_ERROR;
674 sprintf(Finley_ErrorMsg,"Finley does not allow interpolation across domains yet.");
675 checkFinleyError();
676 }
677
678 //
679 // calculates the integral of a function defined of arg:
680 //
681 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const
682 {
683 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
684 if (argDomain!=*this)
685 throw FinleyAdapterException("Error - Illegal domain of integration kernel");
686
687 Finley_Mesh* mesh=m_finleyMesh.get();
688 switch(arg.getFunctionSpace().getTypeCode()) {
689 case(Nodes):
690 Finley_ErrorCode=TYPE_ERROR;
691 sprintf(Finley_ErrorMsg,"Integral of data on nodes is not supported.");
692 break;
693 case(Elements):
694 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);
695 break;
696 case(FaceElements):
697 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);
698 break;
699 case(Points):
700 Finley_ErrorCode=TYPE_ERROR;
701 sprintf(Finley_ErrorMsg,"Integral of data on points is not supported.");
702 break;
703 case(ContactElementsZero):
704 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
705 break;
706 case(ContactElementsOne):
707 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
708 break;
709 case(DegreesOfFreedom):
710 Finley_ErrorCode=TYPE_ERROR;
711 sprintf(Finley_ErrorMsg,"Integral of data on degrees of freedom is not supported.");
712 break;
713 case(ReducedDegreesOfFreedom):
714 Finley_ErrorCode=TYPE_ERROR;
715 sprintf(Finley_ErrorMsg,"Integral of data on reduced degrees of freedom is not supported.");
716 break;
717 default:
718 Finley_ErrorCode=TYPE_ERROR;
719 sprintf(Finley_ErrorMsg,"Integrals: Finley does not know anything about function space type %d",arg.getFunctionSpace().getTypeCode());
720 break;
721 }
722 checkFinleyError();
723 }
724
725 //
726 // calculates the gradient of arg:
727 //
728 void MeshAdapter::setToGradient(Data& grad,const Data& arg) const
729 {
730 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
731 if (argDomain!=*this)
732 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
733 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
734 if (gradDomain!=*this)
735 throw FinleyAdapterException("Error - Illegal domain of gradient");
736
737 Finley_Mesh* mesh=m_finleyMesh.get();
738 switch(grad.getFunctionSpace().getTypeCode()) {
739 case(Nodes):
740 Finley_ErrorCode=TYPE_ERROR;
741 sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");
742 break;
743 case(Elements):
744 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));
745 break;
746 case(FaceElements):
747 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));
748 break;
749 case(Points):
750 Finley_ErrorCode=TYPE_ERROR;
751 sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");
752 break;
753 case(ContactElementsZero):
754 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
755 break;
756 case(ContactElementsOne):
757 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
758 break;
759 case(DegreesOfFreedom):
760 Finley_ErrorCode=TYPE_ERROR;
761 sprintf(Finley_ErrorMsg,"Gradient at degrees of freedom is not supported.");
762 break;
763 case(ReducedDegreesOfFreedom):
764 Finley_ErrorCode=TYPE_ERROR;
765 sprintf(Finley_ErrorMsg,"Gradient at reduced degrees of freedom is not supported.");
766 break;
767 default:
768 Finley_ErrorCode=TYPE_ERROR;
769 sprintf(Finley_ErrorMsg,"Gradient: Finley does not know anything about function space type %d",arg.getFunctionSpace().getTypeCode());
770 break;
771 }
772 checkFinleyError();
773 }
774
775 //
776 // returns the size of elements:
777 //
778 void MeshAdapter::setToSize(Data& size) const
779 {
780 Finley_Mesh* mesh=m_finleyMesh.get();
781 escriptDataC tmp=size.getDataC();
782 switch(size.getFunctionSpace().getTypeCode()) {
783 case(Nodes):
784 Finley_ErrorCode=TYPE_ERROR;
785 sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");
786 break;
787 case(Elements):
788 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
789 break;
790 case(FaceElements):
791 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
792 break;
793 case(Points):
794 Finley_ErrorCode=TYPE_ERROR;
795 sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");
796 break;
797 case(ContactElementsZero):
798 case(ContactElementsOne):
799 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
800 break;
801 case(DegreesOfFreedom):
802 Finley_ErrorCode=TYPE_ERROR;
803 sprintf(Finley_ErrorMsg,"Size of degrees of freedom is not supported.");
804 break;
805 case(ReducedDegreesOfFreedom):
806 Finley_ErrorCode=TYPE_ERROR;
807 sprintf(Finley_ErrorMsg,"Size of reduced degrees of freedom is not supported.");
808 break;
809 default:
810 Finley_ErrorCode=TYPE_ERROR;
811 sprintf(Finley_ErrorMsg,"Element size: Finley does not know anything about function space type %d",size.getFunctionSpace().getTypeCode());
812 break;
813 }
814 checkFinleyError();
815 }
816
817 // sets the location of nodes:
818 void MeshAdapter::setNewX(const Data& new_x)
819 {
820 Finley_Mesh* mesh=m_finleyMesh.get();
821 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
822 if (newDomain!=*this)
823 throw FinleyAdapterException("Error - Illegal domain of new point locations");
824 Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));
825 checkFinleyError();
826 }
827
828 // saves a data array in openDX format:
829 void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const
830 {
831 Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));
832 checkFinleyError();
833 }
834
835 // saves a data array in openVTK format:
836 void MeshAdapter::saveVTK(const std::string& filename,const Data& arg) const
837 {
838 Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));
839 checkFinleyError();
840 }
841
842 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
843 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
844 const int row_blocksize,
845 const escript::FunctionSpace& row_functionspace,
846 const int column_blocksize,
847 const escript::FunctionSpace& column_functionspace,
848 const int type) const
849 {
850 int reduceRowOrder=0;
851 int reduceColOrder=0;
852 // is the domain right?
853 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
854 if (row_domain!=*this)
855 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
856 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
857 if (col_domain!=*this)
858 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
859 // is the function space type right
860 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
861 reduceRowOrder=0;
862 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
863 reduceRowOrder=1;
864 } else {
865 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
866 }
867 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
868 reduceColOrder=0;
869 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
870 reduceColOrder=1;
871 } else {
872 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
873 }
874 // generate matrix:
875
876 Finley_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
877 checkFinleyError();
878 Finley_SystemMatrix* fsystemMatrix=Finley_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
879 checkFinleyError();
880 Finley_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
881 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
882 }
883
884 //
885 // vtkObject MeshAdapter::createVtkObject() const
886 // TODO:
887 //
888
889 //
890 // returns true if data at the atom_type is considered as being cell centered:
891 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
892 {
893 switch(functionSpaceCode) {
894 case(Nodes):
895 case(DegreesOfFreedom):
896 case(ReducedDegreesOfFreedom):
897 return false;
898 break;
899 case(Elements):
900 case(FaceElements):
901 case(Points):
902 case(ContactElementsZero):
903 case(ContactElementsOne):
904 return true;
905 break;
906 default:
907 Finley_ErrorCode=TYPE_ERROR;
908 sprintf(Finley_ErrorMsg,"Cell: Finley does not know anything about function space type %d",functionSpaceCode);
909 break;
910 }
911 checkFinleyError();
912 return false;
913 }
914
915 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
916 {
917 switch(functionSpaceType_source) {
918 case(Nodes):
919 switch(functionSpaceType_target) {
920 case(Nodes):
921 case(ReducedDegreesOfFreedom):
922 case(DegreesOfFreedom):
923 case(Elements):
924 case(FaceElements):
925 case(Points):
926 case(ContactElementsZero):
927 case(ContactElementsOne):
928 return true;
929 default:
930 Finley_ErrorCode=TYPE_ERROR;
931 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);
932 }
933 break;
934 case(Elements):
935 if (functionSpaceType_target==Elements) {
936 return true;
937 } else {
938 return false;
939 }
940 case(FaceElements):
941 if (functionSpaceType_target==FaceElements) {
942 return true;
943 } else {
944 return false;
945 }
946 case(Points):
947 if (functionSpaceType_target==Points) {
948 return true;
949 } else {
950 return false;
951 }
952 case(ContactElementsZero):
953 case(ContactElementsOne):
954 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
955 return true;
956 } else {
957 return false;
958 }
959 case(DegreesOfFreedom):
960 switch(functionSpaceType_target) {
961 case(ReducedDegreesOfFreedom):
962 case(DegreesOfFreedom):
963 case(Nodes):
964 case(Elements):
965 case(FaceElements):
966 case(Points):
967 case(ContactElementsZero):
968 case(ContactElementsOne):
969 return true;
970 default:
971 Finley_ErrorCode=TYPE_ERROR;
972 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);
973 }
974 break;
975 case(ReducedDegreesOfFreedom):
976 switch(functionSpaceType_target) {
977 case(ReducedDegreesOfFreedom):
978 case(Elements):
979 case(FaceElements):
980 case(Points):
981 case(ContactElementsZero):
982 case(ContactElementsOne):
983 return true;
984 case(Nodes):
985 case(DegreesOfFreedom):
986 return false;
987 default:
988 Finley_ErrorCode=TYPE_ERROR;
989 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);
990 }
991 break;
992 default:
993 Finley_ErrorCode=TYPE_ERROR;
994 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_source);
995 break;
996 }
997 checkFinleyError();
998 return false;
999 }
1000
1001 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1002 {
1003 return false;
1004 }
1005
1006 bool MeshAdapter::operator==(const AbstractDomain& other) const
1007 {
1008 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1009 if (temp!=0) {
1010 return (m_finleyMesh==temp->m_finleyMesh);
1011 } else {
1012 return false;
1013 }
1014 }
1015
1016 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1017 {
1018 return !(operator==(other));
1019 }
1020
1021 int MeshAdapter::getSystemMatrixTypeId(const int solver, const bool symmetry) const
1022 {
1023 int out=Finley_SystemMatrix_getSystemMatrixTypeId(solver,symmetry?1:0);
1024 checkFinleyError();
1025 return out;
1026 }
1027
1028 Data MeshAdapter::getX() const
1029 {
1030 return continuousFunction(asAbstractContinuousDomain()).getX();
1031 }
1032
1033 Data MeshAdapter::getNormal() const
1034 {
1035 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1036 }
1037
1038 Data MeshAdapter::getSize() const
1039 {
1040 return function(asAbstractContinuousDomain()).getSize();
1041 }
1042
1043 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1044 {
1045 int* tagList;
1046 int numTags;
1047 getTagList(functionSpaceType, &tagList, &numTags);
1048 return tagList[sampleNo];
1049 }
1050
1051 int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const
1052 {
1053 int* referenceNoList;
1054 int numReferenceNo;
1055 getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);
1056 return referenceNoList[sampleNo];
1057 }
1058
1059 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26