/[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 110 - (hide annotations)
Mon Feb 14 04:14:42 2005 UTC (14 years, 2 months ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 37679 byte(s)
*** empty log message ***

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