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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26