/[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 148 - (hide annotations)
Tue Aug 23 01:24:31 2005 UTC (13 years, 8 months ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 38264 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-08-23

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