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

Annotation of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Tue Oct 26 06:53:54 2004 UTC (14 years, 6 months ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 34463 byte(s)
Initial revision

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26