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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 532 - (hide 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 jgs 102 // $Id$
2 jgs 82 /*
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 jgs 472
16 jgs 203 #include "MeshAdapter.h"
17 jgs 149
18 jgs 480 #include "Data.h"
19     #include "DataFactory.h"
20    
21 jgs 82 using namespace std;
22     using namespace escript;
23    
24     namespace finley {
25 jgs 149
26 jgs 82 //
27 jgs 149 // define the static constants
28 jgs 82 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 jgs 149
47 jgs 82 //
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 jgs 149
134 jgs 82 int MeshAdapter::getFunctionCode() const
135     {
136     return Elements;
137     }
138 jgs 149
139 jgs 82 int MeshAdapter::getFunctionOnBoundaryCode() const
140     {
141     return FaceElements;
142     }
143 jgs 149
144 jgs 82 int MeshAdapter::getFunctionOnContactZeroCode() const
145     {
146     return ContactElementsZero;
147     }
148 jgs 149
149 jgs 82 int MeshAdapter::getFunctionOnContactOneCode() const
150     {
151     return ContactElementsOne;
152     }
153    
154     int MeshAdapter::getSolutionCode() const
155     {
156     return DegreesOfFreedom;
157     }
158 jgs 149
159 jgs 82 int MeshAdapter::getReducedSolutionCode() const
160     {
161     return ReducedDegreesOfFreedom;
162     }
163 jgs 149
164 jgs 82 int MeshAdapter::getDiracDeltaFunctionCode() const
165     {
166     return Points;
167     }
168 jgs 149
169 jgs 82 //
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 jgs 102 *tagList=NULL;
176 jgs 82 *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 jgs 150 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
230 jgs 82 throw FinleyAdapterException(temp.str());
231     break;
232     }
233 jgs 110 if (*tagList==NULL) {
234     stringstream temp;
235     temp << "Error - no tags available for " << functionSpaceType << " for domain: " << getDescription();
236     throw FinleyAdapterException(temp.str());
237     }
238 jgs 82 return;
239     }
240 jgs 149
241 jgs 110 //
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 jgs 150 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
302 jgs 110 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 jgs 82
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 jgs 149
323 jgs 82 //
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 jgs 150 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
382 jgs 82 throw FinleyAdapterException(temp.str());
383     break;
384     }
385     return pair<int,int>(numDataPointsPerSample,numSamples);
386     }
387 jgs 149
388 jgs 82 //
389     // adds linear PDE of second order into a given stiffness matrix and right hand side:
390     //
391     void MeshAdapter::addPDEToSystem(
392 jgs 480 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 jgs 82 {
397     Finley_Mesh* mesh=m_finleyMesh.get();
398 jgs 150 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(),&(rhs.getDataC()),
399 jgs 82 &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));
400     checkFinleyError();
401     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,
402 jgs 150 mat.getPaso_SystemMatrix(),
403 jgs 82 &(rhs.getDataC()),
404     &(d.getDataC()),&(y.getDataC()),
405     Finley_Assemble_handelShapeMissMatch_Mean_out);
406     checkFinleyError();
407 jgs 147 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,
408 jgs 150 mat.getPaso_SystemMatrix(),
409 jgs 82 &(rhs.getDataC()),
410     &(d_contact.getDataC()),
411     &(y_contact.getDataC()),
412     Finley_Assemble_handelShapeMissMatch_Step_out);
413     checkFinleyError();
414     }
415 jgs 149
416 jgs 82 //
417 jgs 102 // adds linear PDE of second order into the right hand side only
418     //
419 jgs 480 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 jgs 102 {
422     Finley_Mesh* mesh=m_finleyMesh.get();
423 jgs 147
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 jgs 102 checkFinleyError();
427 jgs 147
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 jgs 102 checkFinleyError();
432 jgs 147 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 jgs 102 checkFinleyError();
435     }
436 jgs 149
437 jgs 102 //
438 jgs 82 // interpolates data between different function spaces:
439     //
440 jgs 480 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
441 jgs 82 {
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 jgs 150 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 jgs 82 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 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
485 jgs 82 }
486     break;
487     case(FaceElements):
488     if (target.getFunctionSpace().getTypeCode()==FaceElements) {
489     Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
490     } else {
491 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
492 jgs 82 break;
493     }
494     case(Points):
495     if (target.getFunctionSpace().getTypeCode()==Points) {
496     Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
497     } else {
498 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
499 jgs 82 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 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
508 jgs 82 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 jgs 150 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 jgs 82 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 jgs 147 case(Nodes):
557 jgs 150 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
558 jgs 147 break;
559 jgs 82 case(DegreesOfFreedom):
560 jgs 150 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
561 jgs 82 break;
562     default:
563 jgs 150 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 jgs 82 break;
567     }
568     break;
569     default:
570 jgs 150 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 jgs 82 break;
574     }
575     checkFinleyError();
576     }
577    
578     //
579     // copies the locations of sample points into x:
580     //
581 jgs 480 void MeshAdapter::setToX(escript::Data& arg) const
582 jgs 82 {
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 jgs 480 escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
592 jgs 82 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 jgs 149
599 jgs 82 //
600     // return the normal vectors at the location of data points as a Data object:
601     //
602 jgs 480 void MeshAdapter::setToNormal(escript::Data& normal) const
603 jgs 82 {
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 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
611 jgs 82 break;
612     case(Elements):
613 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
614 jgs 82 break;
615     case (FaceElements):
616     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
617     break;
618     case(Points):
619 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
620 jgs 82 break;
621     case (ContactElementsOne):
622     case (ContactElementsZero):
623     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
624     break;
625     case(DegreesOfFreedom):
626 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
627 jgs 82 break;
628     case(ReducedDegreesOfFreedom):
629 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
630 jgs 82 break;
631     default:
632 jgs 150 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 jgs 82 break;
636     }
637     checkFinleyError();
638     }
639 jgs 149
640 jgs 82 //
641     // interpolates data to other domain:
642     //
643 jgs 480 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
644 jgs 82 {
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 jgs 150 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
650 jgs 82 }
651 jgs 149
652 jgs 82 //
653     // calculates the integral of a function defined of arg:
654     //
655 jgs 480 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
656 jgs 82 {
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 jgs 150 throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");
665 jgs 82 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 jgs 150 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
674 jgs 82 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 jgs 150 throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");
683 jgs 82 break;
684     case(ReducedDegreesOfFreedom):
685 jgs 150 throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");
686 jgs 82 break;
687     default:
688 jgs 150 stringstream temp;
689     temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
690     throw FinleyAdapterException(temp.str());
691 jgs 82 break;
692     }
693     checkFinleyError();
694     }
695 jgs 149
696 jgs 82 //
697     // calculates the gradient of arg:
698     //
699 jgs 480 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
700 jgs 82 {
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 jgs 150 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
712 jgs 82 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 jgs 150 throw FinleyAdapterException("Error - Gradient at points is not supported.");
721 jgs 82 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 jgs 150 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
730 jgs 82 break;
731     case(ReducedDegreesOfFreedom):
732 jgs 150 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
733 jgs 82 break;
734     default:
735 jgs 150 stringstream temp;
736     temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
737     throw FinleyAdapterException(temp.str());
738 jgs 82 break;
739     }
740     checkFinleyError();
741     }
742 jgs 149
743 jgs 82 //
744     // returns the size of elements:
745     //
746 jgs 480 void MeshAdapter::setToSize(escript::Data& size) const
747 jgs 82 {
748     Finley_Mesh* mesh=m_finleyMesh.get();
749     escriptDataC tmp=size.getDataC();
750     switch(size.getFunctionSpace().getTypeCode()) {
751     case(Nodes):
752 jgs 150 throw FinleyAdapterException("Error - Size of nodes is not supported.");
753 jgs 82 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 jgs 150 throw FinleyAdapterException("Error - Size of point elements is not supported.");
762 jgs 82 break;
763     case(ContactElementsZero):
764     case(ContactElementsOne):
765     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
766     break;
767     case(DegreesOfFreedom):
768 jgs 150 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
769 jgs 82 break;
770     case(ReducedDegreesOfFreedom):
771 jgs 150 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
772 jgs 82 break;
773     default:
774 jgs 150 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 jgs 82 break;
778     }
779     checkFinleyError();
780     }
781 jgs 149
782 jgs 82 // sets the location of nodes:
783 jgs 480 void MeshAdapter::setNewX(const escript::Data& new_x)
784 jgs 82 {
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 gross 532 Finley_Mesh_setCoordinates(mesh,&(new_x.getDataC()));
790 jgs 82 checkFinleyError();
791     }
792 jgs 149
793 jgs 82 // saves a data array in openDX format:
794 jgs 153 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 jgs 480 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
806 jgs 153 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 jgs 82 }
823 jgs 149
824 jgs 110 // saves a data array in openVTK format:
825 jgs 153 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 jgs 480 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
837 jgs 153 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 jgs 110 }
854 jgs 153
855    
856 jgs 82 // 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 jgs 102 const int type) const
863 jgs 82 {
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 jgs 102
890 jgs 150 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
891 jgs 82 checkFinleyError();
892 jgs 150 Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
893     checkPasoError();
894     Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
895 jgs 82 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
896     }
897 jgs 149
898 jgs 82 //
899     // vtkObject MeshAdapter::createVtkObject() const
900     // TODO:
901     //
902 jgs 149
903 jgs 82 //
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 jgs 150 stringstream temp;
922     temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
923     throw FinleyAdapterException(temp.str());
924 jgs 82 break;
925     }
926     checkFinleyError();
927     return false;
928     }
929 jgs 149
930 jgs 82 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 jgs 150 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 jgs 82 }
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 jgs 150 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 jgs 82 }
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 jgs 147 case(Nodes):
1002 jgs 82 case(DegreesOfFreedom):
1003     return false;
1004     default:
1005 jgs 150 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 jgs 82 }
1009     break;
1010     default:
1011 jgs 150 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 jgs 82 break;
1015     }
1016     checkFinleyError();
1017     return false;
1018     }
1019 jgs 149
1020 jgs 82 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1021     {
1022     return false;
1023     }
1024 jgs 104
1025 jgs 121 bool MeshAdapter::operator==(const AbstractDomain& other) const
1026 jgs 82 {
1027 jgs 121 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 jgs 82 }
1034    
1035 jgs 121 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1036 jgs 104 {
1037 jgs 121 return !(operator==(other));
1038 jgs 104 }
1039    
1040 jgs 150 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1041 jgs 102 {
1042 jgs 150 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1043     checkPasoError();
1044 jgs 102 return out;
1045     }
1046 jgs 149
1047 jgs 480 escript::Data MeshAdapter::getX() const
1048 jgs 102 {
1049     return continuousFunction(asAbstractContinuousDomain()).getX();
1050     }
1051 jgs 149
1052 jgs 480 escript::Data MeshAdapter::getNormal() const
1053 jgs 102 {
1054     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1055     }
1056 jgs 149
1057 jgs 480 escript::Data MeshAdapter::getSize() const
1058 jgs 102 {
1059     return function(asAbstractContinuousDomain()).getSize();
1060     }
1061    
1062 jgs 110 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 jgs 102
1070 jgs 110 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 jgs 82 } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
  ViewVC Help
Powered by ViewVC 1.1.26