/[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 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (14 years, 3 months ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 35297 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 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
184 {
185 //
186 // It is assumed the sampleNo has been checked
187 // before calling this function.
188 int* tagList;
189 int numTags;
190 getTagList(functionSpaceType, &tagList, &numTags);
191 return tagList[sampleNo];
192 }
193
194 //
195 // returns a pointer to the tag list of samples of functionSpaceType
196 //
197 void MeshAdapter::getTagList(int functionSpaceType, int** tagList,
198 int* numTags) const
199 {
200 *tagList=NULL;
201 *numTags=0;
202 Finley_Mesh* mesh=m_finleyMesh.get();
203 switch (functionSpaceType) {
204 case(Nodes):
205 if (mesh->Nodes!=NULL) {
206 *tagList=mesh->Nodes->Tag;
207 *numTags=mesh->Nodes->numNodes;
208 }
209 break;
210 case(Elements):
211 if (mesh->Elements!=NULL) {
212 *tagList=mesh->Elements->Tag;
213 *numTags=mesh->Elements->numElements;
214 }
215 break;
216 case(FaceElements):
217 if (mesh->FaceElements!=NULL) {
218 *tagList=mesh->FaceElements->Tag;
219 *numTags=mesh->FaceElements->numElements;
220 }
221 break;
222 case(Points):
223 if (mesh->Points!=NULL) {
224 *tagList=mesh->Points->Tag;
225 *numTags=mesh->Points->numElements;
226 }
227 break;
228 case(ContactElementsZero):
229 if (mesh->ContactElements!=NULL) {
230 *tagList=mesh->ContactElements->Tag;
231 *numTags=mesh->ContactElements->numElements;
232 }
233 break;
234 case(ContactElementsOne):
235 if (mesh->ContactElements!=NULL) {
236 *tagList=mesh->ContactElements->Tag;
237 *numTags=mesh->ContactElements->numElements;
238 }
239 break;
240 case(DegreesOfFreedom):
241 if (mesh->Nodes!=NULL) {
242 *tagList=NULL;
243 *numTags=0;
244 }
245 break;
246 case(ReducedDegreesOfFreedom):
247 if (mesh->Nodes!=NULL) {
248 *tagList=NULL;
249 *numTags=0;
250 }
251 break;
252 default:
253 stringstream temp;
254 temp << "Error - Invalid function space type: "
255 << functionSpaceType << " for domain: " << getDescription();
256 throw FinleyAdapterException(temp.str());
257 break;
258 }
259 return;
260 }
261
262 //
263 // return the spatial dimension of the Mesh:
264 //
265 int MeshAdapter::getDim() const
266 {
267 int numDim=Finley_Mesh_getDim(m_finleyMesh.get());
268 checkFinleyError();
269 return numDim;
270 }
271 //
272 // return the number of data points per sample and the number of samples
273 // needed to represent data on a parts of the mesh.
274 //
275 pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
276 {
277 int numDataPointsPerSample=0;
278 int numSamples=0;
279 Finley_Mesh* mesh=m_finleyMesh.get();
280 switch (functionSpaceCode) {
281 case(Nodes):
282 numDataPointsPerSample=1;
283 if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
284 break;
285 case(Elements):
286 if (mesh->Elements!=NULL) {
287 numSamples=mesh->Elements->numElements;
288 numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
289 }
290 break;
291 case(FaceElements):
292 if (mesh->FaceElements!=NULL) {
293 numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
294 numSamples=mesh->FaceElements->numElements;
295 }
296 break;
297 case(Points):
298 if (mesh->Points!=NULL) {
299 numDataPointsPerSample=1;
300 numSamples=mesh->Points->numElements;
301 }
302 break;
303 case(ContactElementsZero):
304 if (mesh->ContactElements!=NULL) {
305 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
306 numSamples=mesh->ContactElements->numElements;
307 }
308 break;
309 case(ContactElementsOne):
310 if (mesh->ContactElements!=NULL) {
311 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
312 numSamples=mesh->ContactElements->numElements;
313 }
314 break;
315 case(DegreesOfFreedom):
316 if (mesh->Nodes!=NULL) {
317 numDataPointsPerSample=1;
318 numSamples=mesh->Nodes->numDegreesOfFreedom;
319 }
320 break;
321 case(ReducedDegreesOfFreedom):
322 if (mesh->Nodes!=NULL) {
323 numDataPointsPerSample=1;
324 numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
325 }
326 break;
327 default:
328 stringstream temp;
329 temp << "Error - Invalid function space type: "
330 << functionSpaceCode << " for domain: " << getDescription();
331 throw FinleyAdapterException(temp.str());
332 break;
333 }
334 return pair<int,int>(numDataPointsPerSample,numSamples);
335 }
336 //
337 // adds linear PDE of second order into a given stiffness matrix and right hand side:
338 //
339 void MeshAdapter::addPDEToSystem(
340 SystemMatrixAdapter& mat, Data& rhs,
341 const Data& A, const Data& B, const Data& C,const Data& D,const Data& X,const Data& Y,
342 const Data& d, const Data& y,
343 const Data& d_contact,const Data& y_contact) const
344 {
345 Finley_Mesh* mesh=m_finleyMesh.get();
346 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getFinley_SystemMatrix(),&(rhs.getDataC()),
347 &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));
348 checkFinleyError();
349 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,
350 mat.getFinley_SystemMatrix(),
351 &(rhs.getDataC()),
352 &(d.getDataC()),&(y.getDataC()),
353 Finley_Assemble_handelShapeMissMatch_Mean_out);
354 checkFinleyError();
355 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,
356 mat.getFinley_SystemMatrix(),
357 &(rhs.getDataC()),
358 &(d_contact.getDataC()),
359 &(y_contact.getDataC()),
360 Finley_Assemble_handelShapeMissMatch_Step_out);
361 checkFinleyError();
362 }
363 //
364 // adds linear PDE of second order into the right hand side only
365 //
366 void MeshAdapter::addPDEToRHS( Data& rhs,
367 const Data& X,const Data& Y, const Data& y, const Data& y_contact) const
368 {
369 Finley_Mesh* mesh=m_finleyMesh.get();
370 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));
371 checkFinleyError();
372 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),
373 Finley_Assemble_handelShapeMissMatch_Mean_out);
374 // cout << "Calling :addPDEToRHS." << endl;
375 checkFinleyError();
376 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),
377 Finley_Assemble_handelShapeMissMatch_Step_out);
378 // cout << "Calling :addPDEToRHS." << endl;
379 checkFinleyError();
380 }
381 //
382 // interpolates data between different function spaces:
383 //
384 void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const
385 {
386 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
387 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
388 if (inDomain!=*this)
389 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
390 if (targetDomain!=*this)
391 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
392
393 Finley_Mesh* mesh=m_finleyMesh.get();
394 switch(in.getFunctionSpace().getTypeCode()) {
395 case(Nodes):
396 switch(target.getFunctionSpace().getTypeCode()) {
397 case(Nodes):
398 case(ReducedDegreesOfFreedom):
399 case(DegreesOfFreedom):
400 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
401 break;
402 case(Elements):
403 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
404 break;
405 case(FaceElements):
406 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
407 break;
408 case(Points):
409 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
410 break;
411 case(ContactElementsZero):
412 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
413 break;
414 case(ContactElementsOne):
415 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
416 break;
417 default:
418 Finley_ErrorCode=TYPE_ERROR;
419 sprintf(Finley_ErrorMsg,"Interpolation on Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());
420 break;
421 }
422 break;
423 case(Elements):
424 if (target.getFunctionSpace().getTypeCode()==Elements) {
425 Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));
426 } else {
427 Finley_ErrorCode=TYPE_ERROR;
428 sprintf(Finley_ErrorMsg,"No interpolation with data on elements possible.");
429 }
430 break;
431 case(FaceElements):
432 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
433 Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
434 } else {
435 Finley_ErrorCode=TYPE_ERROR;
436 sprintf(Finley_ErrorMsg,"No interpolation with data on face elements possible.");
437 break;
438 }
439 case(Points):
440 if (target.getFunctionSpace().getTypeCode()==Points) {
441 Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
442 } else {
443 Finley_ErrorCode=TYPE_ERROR;
444 sprintf(Finley_ErrorMsg,"No interpolation with data on points possible.");
445 break;
446 }
447 break;
448 case(ContactElementsZero):
449 case(ContactElementsOne):
450 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
451 Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));
452 } else {
453 Finley_ErrorCode=TYPE_ERROR;
454 sprintf(Finley_ErrorMsg,"No interpolation with data on contact elements possible.");
455 break;
456 }
457 break;
458 case(DegreesOfFreedom):
459 switch(target.getFunctionSpace().getTypeCode()) {
460 case(ReducedDegreesOfFreedom):
461 case(DegreesOfFreedom):
462 case(Nodes):
463 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
464 break;
465 case(Elements):
466 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
467 break;
468 case(FaceElements):
469 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
470 break;
471 case(Points):
472 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
473 break;
474 case(ContactElementsZero):
475 case(ContactElementsOne):
476 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
477 break;
478 default:
479 Finley_ErrorCode=TYPE_ERROR;
480 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());
481 break;
482 }
483 break;
484 case(ReducedDegreesOfFreedom):
485 switch(target.getFunctionSpace().getTypeCode()) {
486 case(ReducedDegreesOfFreedom):
487 case(Nodes):
488 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
489 break;
490 case(Elements):
491 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
492 break;
493 case(FaceElements):
494 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
495 break;
496 case(Points):
497 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
498 break;
499 case(ContactElementsZero):
500 case(ContactElementsOne):
501 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
502 break;
503 case(DegreesOfFreedom):
504 Finley_ErrorCode=TYPE_ERROR;
505 sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
506 break;
507 default:
508 Finley_ErrorCode=TYPE_ERROR;
509 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());
510 break;
511 }
512 break;
513 default:
514 Finley_ErrorCode=TYPE_ERROR;
515 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",in.getFunctionSpace().getTypeCode());
516 break;
517 }
518 checkFinleyError();
519 }
520
521 //
522 // copies the locations of sample points into x:
523 //
524 void MeshAdapter::setToX(Data& arg) const
525 {
526 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
527 if (argDomain!=*this)
528 throw FinleyAdapterException("Error - Illegal domain of data point locations");
529 Finley_Mesh* mesh=m_finleyMesh.get();
530 // in case of values node coordinates we can do the job directly:
531 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
532 Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
533 } else {
534 Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
535 Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
536 // this is then interpolated onto arg:
537 interpolateOnDomain(arg,tmp_data);
538 }
539 checkFinleyError();
540 }
541 //
542 // return the normal vectors at the location of data points as a Data object:
543 //
544 void MeshAdapter::setToNormal(Data& normal) const
545 {
546 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
547 if (normalDomain!=*this)
548 throw FinleyAdapterException("Error - Illegal domain of normal locations");
549 Finley_Mesh* mesh=m_finleyMesh.get();
550 switch(normal.getFunctionSpace().getTypeCode()) {
551 case(Nodes):
552 Finley_ErrorCode=TYPE_ERROR;
553 sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for nodes");
554 break;
555 case(Elements):
556 Finley_ErrorCode=TYPE_ERROR;
557 sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for elements");
558 break;
559 case (FaceElements):
560 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
561 break;
562 case(Points):
563 Finley_ErrorCode=TYPE_ERROR;
564 sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for point elements");
565 break;
566 case (ContactElementsOne):
567 case (ContactElementsZero):
568 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
569 break;
570 case(DegreesOfFreedom):
571 Finley_ErrorCode=TYPE_ERROR;
572 sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for degrees of freedom.");
573 break;
574 case(ReducedDegreesOfFreedom):
575 Finley_ErrorCode=TYPE_ERROR;
576 sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for reduced degrees of freedom.");
577 break;
578 default:
579 Finley_ErrorCode=TYPE_ERROR;
580 sprintf(Finley_ErrorMsg,"Normal Vectors: Finley does not know anything about function space type %d",normal.getFunctionSpace().getTypeCode());
581 break;
582 }
583 checkFinleyError();
584 }
585 //
586 // interpolates data to other domain:
587 //
588 void MeshAdapter::interpolateACross(Data& target,const Data& source) const
589 {
590 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
591 if (targetDomain!=*this)
592 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
593
594 Finley_ErrorCode=SYSTEM_ERROR;
595 sprintf(Finley_ErrorMsg,"Finley does not allow interpolation across domains yet.");
596 checkFinleyError();
597 }
598 //
599 // calculates the integral of a function defined of arg:
600 //
601 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const 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 integration kernel");
606
607 Finley_Mesh* mesh=m_finleyMesh.get();
608 switch(arg.getFunctionSpace().getTypeCode()) {
609 case(Nodes):
610 Finley_ErrorCode=TYPE_ERROR;
611 sprintf(Finley_ErrorMsg,"Integral of data on nodes is not supported.");
612 break;
613 case(Elements):
614 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);
615 break;
616 case(FaceElements):
617 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);
618 break;
619 case(Points):
620 Finley_ErrorCode=TYPE_ERROR;
621 sprintf(Finley_ErrorMsg,"Integral of data on points is not supported.");
622 break;
623 case(ContactElementsZero):
624 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
625 break;
626 case(ContactElementsOne):
627 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
628 break;
629 case(DegreesOfFreedom):
630 Finley_ErrorCode=TYPE_ERROR;
631 sprintf(Finley_ErrorMsg,"Integral of data on degrees of freedom is not supported.");
632 break;
633 case(ReducedDegreesOfFreedom):
634 Finley_ErrorCode=TYPE_ERROR;
635 sprintf(Finley_ErrorMsg,"Integral of data on reduced degrees of freedom is not supported.");
636 break;
637 default:
638 Finley_ErrorCode=TYPE_ERROR;
639 sprintf(Finley_ErrorMsg,"Integrals: Finley does not know anything about function space type %d",arg.getFunctionSpace().getTypeCode());
640 break;
641 }
642 checkFinleyError();
643 }
644 //
645 // calculates the gradient of arg:
646 //
647 void MeshAdapter::setToGradient(Data& grad,const Data& arg) const
648 {
649 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
650 if (argDomain!=*this)
651 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
652 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
653 if (gradDomain!=*this)
654 throw FinleyAdapterException("Error - Illegal domain of gradient");
655
656 Finley_Mesh* mesh=m_finleyMesh.get();
657 switch(grad.getFunctionSpace().getTypeCode()) {
658 case(Nodes):
659 Finley_ErrorCode=TYPE_ERROR;
660 sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");
661 break;
662 case(Elements):
663 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));
664 break;
665 case(FaceElements):
666 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));
667 break;
668 case(Points):
669 Finley_ErrorCode=TYPE_ERROR;
670 sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");
671 break;
672 case(ContactElementsZero):
673 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
674 break;
675 case(ContactElementsOne):
676 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
677 break;
678 case(DegreesOfFreedom):
679 Finley_ErrorCode=TYPE_ERROR;
680 sprintf(Finley_ErrorMsg,"Gradient at degrees of freedom is not supported.");
681 break;
682 case(ReducedDegreesOfFreedom):
683 Finley_ErrorCode=TYPE_ERROR;
684 sprintf(Finley_ErrorMsg,"Gradient at reduced degrees of freedom is not supported.");
685 break;
686 default:
687 Finley_ErrorCode=TYPE_ERROR;
688 sprintf(Finley_ErrorMsg,"Gradient: Finley does not know anything about function space type %d",arg.getFunctionSpace().getTypeCode());
689 break;
690 }
691 checkFinleyError();
692 }
693 //
694 // returns the size of elements:
695 //
696 void MeshAdapter::setToSize(Data& size) const
697 {
698 Finley_Mesh* mesh=m_finleyMesh.get();
699 escriptDataC tmp=size.getDataC();
700 switch(size.getFunctionSpace().getTypeCode()) {
701 case(Nodes):
702 Finley_ErrorCode=TYPE_ERROR;
703 sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");
704 break;
705 case(Elements):
706 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
707 break;
708 case(FaceElements):
709 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
710 break;
711 case(Points):
712 Finley_ErrorCode=TYPE_ERROR;
713 sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");
714 break;
715 case(ContactElementsZero):
716 case(ContactElementsOne):
717 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
718 break;
719 case(DegreesOfFreedom):
720 Finley_ErrorCode=TYPE_ERROR;
721 sprintf(Finley_ErrorMsg,"Size of degrees of freedom is not supported.");
722 break;
723 case(ReducedDegreesOfFreedom):
724 Finley_ErrorCode=TYPE_ERROR;
725 sprintf(Finley_ErrorMsg,"Size of reduced degrees of freedom is not supported.");
726 break;
727 default:
728 Finley_ErrorCode=TYPE_ERROR;
729 sprintf(Finley_ErrorMsg,"Element size: Finley does not know anything about function space type %d",size.getFunctionSpace().getTypeCode());
730 break;
731 }
732 checkFinleyError();
733 }
734 // sets the location of nodes:
735 void MeshAdapter::setNewX(const Data& new_x)
736 {
737 Finley_Mesh* mesh=m_finleyMesh.get();
738 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
739 if (newDomain!=*this)
740 throw FinleyAdapterException("Error - Illegal domain of new point locations");
741 Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));
742 checkFinleyError();
743 }
744 // saves a data array in openDX format:
745 void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const
746 {
747 Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));
748 checkFinleyError();
749 }
750 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
751 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
752 const int row_blocksize,
753 const escript::FunctionSpace& row_functionspace,
754 const int column_blocksize,
755 const escript::FunctionSpace& column_functionspace,
756 const int type) const
757 {
758 int reduceRowOrder=0;
759 int reduceColOrder=0;
760 // is the domain right?
761 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
762 if (row_domain!=*this)
763 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
764 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
765 if (col_domain!=*this)
766 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
767 // is the function space type right
768 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
769 reduceRowOrder=0;
770 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
771 reduceRowOrder=1;
772 } else {
773 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
774 }
775 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
776 reduceColOrder=0;
777 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
778 reduceColOrder=1;
779 } else {
780 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
781 }
782 // generate matrix:
783
784 Finley_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
785 checkFinleyError();
786 Finley_SystemMatrix* fsystemMatrix=Finley_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
787 checkFinleyError();
788 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
789 }
790 //
791 // vtkObject MeshAdapter::createVtkObject() const
792 // TODO:
793 //
794 //
795 // returns true if data at the atom_type is considered as being cell centered:
796 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
797 {
798 switch(functionSpaceCode) {
799 case(Nodes):
800 case(DegreesOfFreedom):
801 case(ReducedDegreesOfFreedom):
802 return false;
803 break;
804 case(Elements):
805 case(FaceElements):
806 case(Points):
807 case(ContactElementsZero):
808 case(ContactElementsOne):
809 return true;
810 break;
811 default:
812 Finley_ErrorCode=TYPE_ERROR;
813 sprintf(Finley_ErrorMsg,"Cell: Finley does not know anything about function space type %d",functionSpaceCode);
814 break;
815 }
816 checkFinleyError();
817 return false;
818 }
819 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
820 {
821 switch(functionSpaceType_source) {
822 case(Nodes):
823 switch(functionSpaceType_target) {
824 case(Nodes):
825 case(ReducedDegreesOfFreedom):
826 case(DegreesOfFreedom):
827 case(Elements):
828 case(FaceElements):
829 case(Points):
830 case(ContactElementsZero):
831 case(ContactElementsOne):
832 return true;
833 default:
834 Finley_ErrorCode=TYPE_ERROR;
835 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);
836 }
837 break;
838 case(Elements):
839 if (functionSpaceType_target==Elements) {
840 return true;
841 } else {
842 return false;
843 }
844 case(FaceElements):
845 if (functionSpaceType_target==FaceElements) {
846 return true;
847 } else {
848 return false;
849 }
850 case(Points):
851 if (functionSpaceType_target==Points) {
852 return true;
853 } else {
854 return false;
855 }
856 case(ContactElementsZero):
857 case(ContactElementsOne):
858 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
859 return true;
860 } else {
861 return false;
862 }
863 case(DegreesOfFreedom):
864 switch(functionSpaceType_target) {
865 case(ReducedDegreesOfFreedom):
866 case(DegreesOfFreedom):
867 case(Nodes):
868 case(Elements):
869 case(FaceElements):
870 case(Points):
871 case(ContactElementsZero):
872 case(ContactElementsOne):
873 return true;
874 default:
875 Finley_ErrorCode=TYPE_ERROR;
876 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);
877 }
878 break;
879 case(ReducedDegreesOfFreedom):
880 switch(functionSpaceType_target) {
881 case(ReducedDegreesOfFreedom):
882 case(Nodes):
883 case(Elements):
884 case(FaceElements):
885 case(Points):
886 case(ContactElementsZero):
887 case(ContactElementsOne):
888 return true;
889 case(DegreesOfFreedom):
890 return false;
891 default:
892 Finley_ErrorCode=TYPE_ERROR;
893 sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);
894 }
895 break;
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_source);
899 break;
900 }
901 checkFinleyError();
902 return false;
903 }
904 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
905 {
906 return false;
907 }
908 bool MeshAdapter::operator==(const MeshAdapter& other) const
909 {
910 return (m_finleyMesh==other.m_finleyMesh);
911 }
912
913 int MeshAdapter::getSystemMatrixTypeId(const int solver, const bool symmetry) const
914 {
915 int out=Finley_SystemMatrix_getSystemMatrixTypeId(solver,symmetry?1:0);
916 checkFinleyError();
917 return out;
918 }
919 Data MeshAdapter::getX() const
920 {
921 return continuousFunction(asAbstractContinuousDomain()).getX();
922 }
923 Data MeshAdapter::getNormal() const
924 {
925 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
926 }
927 Data MeshAdapter::getSize() const
928 {
929 return function(asAbstractContinuousDomain()).getSize();
930 }
931
932
933 bool MeshAdapter::operator!=(const MeshAdapter& other) const
934 {
935 return !operator==(other);
936 }
937 // bool MeshAdapter::operator==(const AbstractDomain& other) const
938 // {
939 // try {
940 // const MeshAdapter& ref = dynamic_cast<const MeshAdapter&>(other);
941 // return (operator==(ref));
942 // }
943 // catch (bad_cast) {
944 // return false;
945 // }
946 // }
947 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26