/[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 472 - (hide annotations)
Fri Jan 27 01:50:59 2006 UTC (13 years, 4 months ago) by jgs
File size: 39959 byte(s)
rationalise all #includes

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