/[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 123 - (show annotations)
Fri Jul 8 04:08:13 2005 UTC (14 years, 2 months ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 37731 byte(s)
Merge of development branch back to main trunk on 2005-07-08

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26