/[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 110 - (show annotations)
Mon Feb 14 04:14:42 2005 UTC (14 years, 1 month ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 37679 byte(s)
*** empty log message ***

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26