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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 100 - (show annotations)
Wed Dec 15 03:48:48 2004 UTC (14 years, 11 months ago) by jgs
File size: 34463 byte(s)
*** empty log message ***

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26