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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 102 - (hide annotations)
Wed Dec 15 07:08:39 2004 UTC (14 years, 4 months ago) by jgs
File size: 35297 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     int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
184     {
185     //
186     // It is assumed the sampleNo has been checked
187     // before calling this function.
188     int* tagList;
189     int numTags;
190     getTagList(functionSpaceType, &tagList, &numTags);
191     return tagList[sampleNo];
192     }
193    
194     //
195     // returns a pointer to the tag list of samples of functionSpaceType
196     //
197     void MeshAdapter::getTagList(int functionSpaceType, int** tagList,
198     int* numTags) const
199     {
200 jgs 102 *tagList=NULL;
201 jgs 82 *numTags=0;
202     Finley_Mesh* mesh=m_finleyMesh.get();
203     switch (functionSpaceType) {
204     case(Nodes):
205     if (mesh->Nodes!=NULL) {
206     *tagList=mesh->Nodes->Tag;
207     *numTags=mesh->Nodes->numNodes;
208     }
209     break;
210     case(Elements):
211     if (mesh->Elements!=NULL) {
212     *tagList=mesh->Elements->Tag;
213     *numTags=mesh->Elements->numElements;
214     }
215     break;
216     case(FaceElements):
217     if (mesh->FaceElements!=NULL) {
218     *tagList=mesh->FaceElements->Tag;
219     *numTags=mesh->FaceElements->numElements;
220     }
221     break;
222     case(Points):
223     if (mesh->Points!=NULL) {
224     *tagList=mesh->Points->Tag;
225     *numTags=mesh->Points->numElements;
226     }
227     break;
228     case(ContactElementsZero):
229     if (mesh->ContactElements!=NULL) {
230     *tagList=mesh->ContactElements->Tag;
231     *numTags=mesh->ContactElements->numElements;
232     }
233     break;
234     case(ContactElementsOne):
235     if (mesh->ContactElements!=NULL) {
236     *tagList=mesh->ContactElements->Tag;
237     *numTags=mesh->ContactElements->numElements;
238     }
239     break;
240     case(DegreesOfFreedom):
241     if (mesh->Nodes!=NULL) {
242     *tagList=NULL;
243     *numTags=0;
244     }
245     break;
246     case(ReducedDegreesOfFreedom):
247     if (mesh->Nodes!=NULL) {
248     *tagList=NULL;
249     *numTags=0;
250     }
251     break;
252     default:
253     stringstream temp;
254     temp << "Error - Invalid function space type: "
255     << functionSpaceType << " for domain: " << getDescription();
256     throw FinleyAdapterException(temp.str());
257     break;
258     }
259     return;
260     }
261    
262     //
263     // return the spatial dimension of the Mesh:
264     //
265     int MeshAdapter::getDim() const
266     {
267     int numDim=Finley_Mesh_getDim(m_finleyMesh.get());
268     checkFinleyError();
269     return numDim;
270     }
271     //
272     // return the number of data points per sample and the number of samples
273     // needed to represent data on a parts of the mesh.
274     //
275     pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
276     {
277     int numDataPointsPerSample=0;
278     int numSamples=0;
279     Finley_Mesh* mesh=m_finleyMesh.get();
280     switch (functionSpaceCode) {
281     case(Nodes):
282     numDataPointsPerSample=1;
283     if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
284     break;
285     case(Elements):
286     if (mesh->Elements!=NULL) {
287     numSamples=mesh->Elements->numElements;
288     numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
289     }
290     break;
291     case(FaceElements):
292     if (mesh->FaceElements!=NULL) {
293     numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
294     numSamples=mesh->FaceElements->numElements;
295     }
296     break;
297     case(Points):
298     if (mesh->Points!=NULL) {
299     numDataPointsPerSample=1;
300     numSamples=mesh->Points->numElements;
301     }
302     break;
303     case(ContactElementsZero):
304     if (mesh->ContactElements!=NULL) {
305     numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
306     numSamples=mesh->ContactElements->numElements;
307     }
308     break;
309     case(ContactElementsOne):
310     if (mesh->ContactElements!=NULL) {
311     numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
312     numSamples=mesh->ContactElements->numElements;
313     }
314     break;
315     case(DegreesOfFreedom):
316     if (mesh->Nodes!=NULL) {
317     numDataPointsPerSample=1;
318     numSamples=mesh->Nodes->numDegreesOfFreedom;
319     }
320     break;
321     case(ReducedDegreesOfFreedom):
322     if (mesh->Nodes!=NULL) {
323     numDataPointsPerSample=1;
324     numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
325     }
326     break;
327     default:
328     stringstream temp;
329     temp << "Error - Invalid function space type: "
330     << functionSpaceCode << " for domain: " << getDescription();
331     throw FinleyAdapterException(temp.str());
332     break;
333     }
334     return pair<int,int>(numDataPointsPerSample,numSamples);
335     }
336     //
337     // adds linear PDE of second order into a given stiffness matrix and right hand side:
338     //
339     void MeshAdapter::addPDEToSystem(
340     SystemMatrixAdapter& mat, Data& rhs,
341 jgs 102 const Data& A, const Data& B, const Data& C,const Data& D,const Data& X,const Data& Y,
342     const Data& d, const Data& y,
343     const Data& d_contact,const Data& y_contact) const
344 jgs 82 {
345     Finley_Mesh* mesh=m_finleyMesh.get();
346     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getFinley_SystemMatrix(),&(rhs.getDataC()),
347     &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));
348     checkFinleyError();
349     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,
350     mat.getFinley_SystemMatrix(),
351     &(rhs.getDataC()),
352     &(d.getDataC()),&(y.getDataC()),
353     Finley_Assemble_handelShapeMissMatch_Mean_out);
354     checkFinleyError();
355     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,
356     mat.getFinley_SystemMatrix(),
357     &(rhs.getDataC()),
358     &(d_contact.getDataC()),
359     &(y_contact.getDataC()),
360     Finley_Assemble_handelShapeMissMatch_Step_out);
361     checkFinleyError();
362     }
363     //
364 jgs 102 // adds linear PDE of second order into the right hand side only
365     //
366     void MeshAdapter::addPDEToRHS( Data& rhs,
367     const Data& X,const Data& Y, const Data& y, const Data& y_contact) const
368     {
369     Finley_Mesh* mesh=m_finleyMesh.get();
370     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));
371     checkFinleyError();
372     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),
373     Finley_Assemble_handelShapeMissMatch_Mean_out);
374     // cout << "Calling :addPDEToRHS." << endl;
375     checkFinleyError();
376     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),
377     Finley_Assemble_handelShapeMissMatch_Step_out);
378     // cout << "Calling :addPDEToRHS." << endl;
379     checkFinleyError();
380     }
381     //
382 jgs 82 // interpolates data between different function spaces:
383     //
384     void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const
385     {
386     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
387     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
388     if (inDomain!=*this)
389     throw FinleyAdapterException("Error - Illegal domain of interpolant.");
390     if (targetDomain!=*this)
391     throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
392    
393     Finley_Mesh* mesh=m_finleyMesh.get();
394     switch(in.getFunctionSpace().getTypeCode()) {
395     case(Nodes):
396     switch(target.getFunctionSpace().getTypeCode()) {
397     case(Nodes):
398     case(ReducedDegreesOfFreedom):
399     case(DegreesOfFreedom):
400     Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
401     break;
402     case(Elements):
403     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
404     break;
405     case(FaceElements):
406     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
407     break;
408     case(Points):
409     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
410     break;
411     case(ContactElementsZero):
412     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
413     break;
414     case(ContactElementsOne):
415     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
416     break;
417     default:
418     Finley_ErrorCode=TYPE_ERROR;
419     sprintf(Finley_ErrorMsg,"Interpolation on Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());
420     break;
421     }
422     break;
423     case(Elements):
424     if (target.getFunctionSpace().getTypeCode()==Elements) {
425     Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));
426     } else {
427     Finley_ErrorCode=TYPE_ERROR;
428     sprintf(Finley_ErrorMsg,"No interpolation with data on elements possible.");
429     }
430     break;
431     case(FaceElements):
432     if (target.getFunctionSpace().getTypeCode()==FaceElements) {
433     Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
434     } else {
435     Finley_ErrorCode=TYPE_ERROR;
436     sprintf(Finley_ErrorMsg,"No interpolation with data on face elements possible.");
437     break;
438     }
439     case(Points):
440     if (target.getFunctionSpace().getTypeCode()==Points) {
441     Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
442     } else {
443     Finley_ErrorCode=TYPE_ERROR;
444     sprintf(Finley_ErrorMsg,"No interpolation with data on points possible.");
445     break;
446     }
447     break;
448     case(ContactElementsZero):
449     case(ContactElementsOne):
450     if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
451     Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));
452     } else {
453     Finley_ErrorCode=TYPE_ERROR;
454     sprintf(Finley_ErrorMsg,"No interpolation with data on contact elements possible.");
455     break;
456     }
457     break;
458     case(DegreesOfFreedom):
459     switch(target.getFunctionSpace().getTypeCode()) {
460     case(ReducedDegreesOfFreedom):
461     case(DegreesOfFreedom):
462     case(Nodes):
463     Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
464     break;
465     case(Elements):
466     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
467     break;
468     case(FaceElements):
469     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
470     break;
471     case(Points):
472     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
473     break;
474     case(ContactElementsZero):
475     case(ContactElementsOne):
476     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
477     break;
478     default:
479     Finley_ErrorCode=TYPE_ERROR;
480     sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());
481     break;
482     }
483     break;
484     case(ReducedDegreesOfFreedom):
485     switch(target.getFunctionSpace().getTypeCode()) {
486     case(ReducedDegreesOfFreedom):
487     case(Nodes):
488     Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
489     break;
490     case(Elements):
491     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
492     break;
493     case(FaceElements):
494     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
495     break;
496     case(Points):
497     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
498     break;
499     case(ContactElementsZero):
500     case(ContactElementsOne):
501     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
502     break;
503     case(DegreesOfFreedom):
504     Finley_ErrorCode=TYPE_ERROR;
505     sprintf(Finley_ErrorMsg,"Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
506     break;
507     default:
508     Finley_ErrorCode=TYPE_ERROR;
509     sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",target.getFunctionSpace().getTypeCode());
510     break;
511     }
512     break;
513     default:
514     Finley_ErrorCode=TYPE_ERROR;
515     sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",in.getFunctionSpace().getTypeCode());
516     break;
517     }
518     checkFinleyError();
519     }
520    
521     //
522     // copies the locations of sample points into x:
523     //
524     void MeshAdapter::setToX(Data& arg) const
525     {
526     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
527     if (argDomain!=*this)
528     throw FinleyAdapterException("Error - Illegal domain of data point locations");
529     Finley_Mesh* mesh=m_finleyMesh.get();
530     // in case of values node coordinates we can do the job directly:
531     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
532     Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
533     } else {
534     Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
535     Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
536     // this is then interpolated onto arg:
537     interpolateOnDomain(arg,tmp_data);
538     }
539     checkFinleyError();
540     }
541     //
542     // return the normal vectors at the location of data points as a Data object:
543     //
544     void MeshAdapter::setToNormal(Data& normal) const
545     {
546     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
547     if (normalDomain!=*this)
548     throw FinleyAdapterException("Error - Illegal domain of normal locations");
549     Finley_Mesh* mesh=m_finleyMesh.get();
550     switch(normal.getFunctionSpace().getTypeCode()) {
551     case(Nodes):
552     Finley_ErrorCode=TYPE_ERROR;
553     sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for nodes");
554     break;
555     case(Elements):
556     Finley_ErrorCode=TYPE_ERROR;
557     sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for elements");
558     break;
559     case (FaceElements):
560     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
561     break;
562     case(Points):
563     Finley_ErrorCode=TYPE_ERROR;
564     sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for point elements");
565     break;
566     case (ContactElementsOne):
567     case (ContactElementsZero):
568     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
569     break;
570     case(DegreesOfFreedom):
571     Finley_ErrorCode=TYPE_ERROR;
572     sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for degrees of freedom.");
573     break;
574     case(ReducedDegreesOfFreedom):
575     Finley_ErrorCode=TYPE_ERROR;
576     sprintf(Finley_ErrorMsg,"Finley does not support surface normal vectors for reduced degrees of freedom.");
577     break;
578     default:
579     Finley_ErrorCode=TYPE_ERROR;
580     sprintf(Finley_ErrorMsg,"Normal Vectors: Finley does not know anything about function space type %d",normal.getFunctionSpace().getTypeCode());
581     break;
582     }
583     checkFinleyError();
584     }
585     //
586     // interpolates data to other domain:
587     //
588     void MeshAdapter::interpolateACross(Data& target,const Data& source) const
589     {
590     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
591     if (targetDomain!=*this)
592     throw FinleyAdapterException("Error - Illegal domain of interpolation target");
593    
594     Finley_ErrorCode=SYSTEM_ERROR;
595     sprintf(Finley_ErrorMsg,"Finley does not allow interpolation across domains yet.");
596     checkFinleyError();
597     }
598     //
599     // calculates the integral of a function defined of arg:
600     //
601     void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const
602     {
603     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
604     if (argDomain!=*this)
605     throw FinleyAdapterException("Error - Illegal domain of integration kernel");
606    
607     Finley_Mesh* mesh=m_finleyMesh.get();
608     switch(arg.getFunctionSpace().getTypeCode()) {
609     case(Nodes):
610     Finley_ErrorCode=TYPE_ERROR;
611     sprintf(Finley_ErrorMsg,"Integral of data on nodes is not supported.");
612     break;
613     case(Elements):
614     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);
615     break;
616     case(FaceElements):
617     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);
618     break;
619     case(Points):
620     Finley_ErrorCode=TYPE_ERROR;
621     sprintf(Finley_ErrorMsg,"Integral of data on points is not supported.");
622     break;
623     case(ContactElementsZero):
624     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
625     break;
626     case(ContactElementsOne):
627     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
628     break;
629     case(DegreesOfFreedom):
630     Finley_ErrorCode=TYPE_ERROR;
631     sprintf(Finley_ErrorMsg,"Integral of data on degrees of freedom is not supported.");
632     break;
633     case(ReducedDegreesOfFreedom):
634     Finley_ErrorCode=TYPE_ERROR;
635     sprintf(Finley_ErrorMsg,"Integral of data on reduced degrees of freedom is not supported.");
636     break;
637     default:
638     Finley_ErrorCode=TYPE_ERROR;
639     sprintf(Finley_ErrorMsg,"Integrals: Finley does not know anything about function space type %d",arg.getFunctionSpace().getTypeCode());
640     break;
641     }
642     checkFinleyError();
643     }
644     //
645     // calculates the gradient of arg:
646     //
647     void MeshAdapter::setToGradient(Data& grad,const Data& arg) const
648     {
649     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
650     if (argDomain!=*this)
651     throw FinleyAdapterException("Error - Illegal domain of gradient argument");
652     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
653     if (gradDomain!=*this)
654     throw FinleyAdapterException("Error - Illegal domain of gradient");
655    
656     Finley_Mesh* mesh=m_finleyMesh.get();
657     switch(grad.getFunctionSpace().getTypeCode()) {
658     case(Nodes):
659     Finley_ErrorCode=TYPE_ERROR;
660     sprintf(Finley_ErrorMsg,"Gradient at nodes is not supported.");
661     break;
662     case(Elements):
663     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));
664     break;
665     case(FaceElements):
666     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));
667     break;
668     case(Points):
669     Finley_ErrorCode=TYPE_ERROR;
670     sprintf(Finley_ErrorMsg,"Gradient at points is not supported.");
671     break;
672     case(ContactElementsZero):
673     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
674     break;
675     case(ContactElementsOne):
676     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
677     break;
678     case(DegreesOfFreedom):
679     Finley_ErrorCode=TYPE_ERROR;
680     sprintf(Finley_ErrorMsg,"Gradient at degrees of freedom is not supported.");
681     break;
682     case(ReducedDegreesOfFreedom):
683     Finley_ErrorCode=TYPE_ERROR;
684     sprintf(Finley_ErrorMsg,"Gradient at reduced degrees of freedom is not supported.");
685     break;
686     default:
687     Finley_ErrorCode=TYPE_ERROR;
688     sprintf(Finley_ErrorMsg,"Gradient: Finley does not know anything about function space type %d",arg.getFunctionSpace().getTypeCode());
689     break;
690     }
691     checkFinleyError();
692     }
693     //
694     // returns the size of elements:
695     //
696     void MeshAdapter::setToSize(Data& size) const
697     {
698     Finley_Mesh* mesh=m_finleyMesh.get();
699     escriptDataC tmp=size.getDataC();
700     switch(size.getFunctionSpace().getTypeCode()) {
701     case(Nodes):
702     Finley_ErrorCode=TYPE_ERROR;
703     sprintf(Finley_ErrorMsg,"Size of nodes is not supported.");
704     break;
705     case(Elements):
706     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
707     break;
708     case(FaceElements):
709     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
710     break;
711     case(Points):
712     Finley_ErrorCode=TYPE_ERROR;
713     sprintf(Finley_ErrorMsg,"Size of point elements is not supported.");
714     break;
715     case(ContactElementsZero):
716     case(ContactElementsOne):
717     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
718     break;
719     case(DegreesOfFreedom):
720     Finley_ErrorCode=TYPE_ERROR;
721     sprintf(Finley_ErrorMsg,"Size of degrees of freedom is not supported.");
722     break;
723     case(ReducedDegreesOfFreedom):
724     Finley_ErrorCode=TYPE_ERROR;
725     sprintf(Finley_ErrorMsg,"Size of reduced degrees of freedom is not supported.");
726     break;
727     default:
728     Finley_ErrorCode=TYPE_ERROR;
729     sprintf(Finley_ErrorMsg,"Element size: Finley does not know anything about function space type %d",size.getFunctionSpace().getTypeCode());
730     break;
731     }
732     checkFinleyError();
733     }
734     // sets the location of nodes:
735     void MeshAdapter::setNewX(const Data& new_x)
736     {
737     Finley_Mesh* mesh=m_finleyMesh.get();
738     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
739     if (newDomain!=*this)
740     throw FinleyAdapterException("Error - Illegal domain of new point locations");
741     Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));
742     checkFinleyError();
743     }
744     // saves a data array in openDX format:
745     void MeshAdapter::saveDX(const std::string& filename,const Data& arg) const
746     {
747     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),&(arg.getDataC()));
748     checkFinleyError();
749     }
750     // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
751     SystemMatrixAdapter MeshAdapter::newSystemMatrix(
752     const int row_blocksize,
753     const escript::FunctionSpace& row_functionspace,
754     const int column_blocksize,
755     const escript::FunctionSpace& column_functionspace,
756 jgs 102 const int type) const
757 jgs 82 {
758     int reduceRowOrder=0;
759     int reduceColOrder=0;
760     // is the domain right?
761     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
762     if (row_domain!=*this)
763     throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
764     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
765     if (col_domain!=*this)
766     throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
767     // is the function space type right
768     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
769     reduceRowOrder=0;
770     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
771     reduceRowOrder=1;
772     } else {
773     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
774     }
775     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
776     reduceColOrder=0;
777     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
778     reduceColOrder=1;
779     } else {
780     throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
781     }
782     // generate matrix:
783 jgs 102
784     Finley_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
785 jgs 82 checkFinleyError();
786 jgs 102 Finley_SystemMatrix* fsystemMatrix=Finley_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
787     checkFinleyError();
788 jgs 82 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
789     }
790     //
791     // vtkObject MeshAdapter::createVtkObject() const
792     // TODO:
793     //
794     //
795     // returns true if data at the atom_type is considered as being cell centered:
796     bool MeshAdapter::isCellOriented(int functionSpaceCode) const
797     {
798     switch(functionSpaceCode) {
799     case(Nodes):
800     case(DegreesOfFreedom):
801     case(ReducedDegreesOfFreedom):
802     return false;
803     break;
804     case(Elements):
805     case(FaceElements):
806     case(Points):
807     case(ContactElementsZero):
808     case(ContactElementsOne):
809     return true;
810     break;
811     default:
812     Finley_ErrorCode=TYPE_ERROR;
813     sprintf(Finley_ErrorMsg,"Cell: Finley does not know anything about function space type %d",functionSpaceCode);
814     break;
815     }
816     checkFinleyError();
817     return false;
818     }
819     bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
820     {
821     switch(functionSpaceType_source) {
822     case(Nodes):
823     switch(functionSpaceType_target) {
824     case(Nodes):
825     case(ReducedDegreesOfFreedom):
826     case(DegreesOfFreedom):
827     case(Elements):
828     case(FaceElements):
829     case(Points):
830     case(ContactElementsZero):
831     case(ContactElementsOne):
832     return true;
833     default:
834     Finley_ErrorCode=TYPE_ERROR;
835     sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);
836     }
837     break;
838     case(Elements):
839     if (functionSpaceType_target==Elements) {
840     return true;
841     } else {
842     return false;
843     }
844     case(FaceElements):
845     if (functionSpaceType_target==FaceElements) {
846     return true;
847     } else {
848     return false;
849     }
850     case(Points):
851     if (functionSpaceType_target==Points) {
852     return true;
853     } else {
854     return false;
855     }
856     case(ContactElementsZero):
857     case(ContactElementsOne):
858     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
859     return true;
860     } else {
861     return false;
862     }
863     case(DegreesOfFreedom):
864     switch(functionSpaceType_target) {
865     case(ReducedDegreesOfFreedom):
866     case(DegreesOfFreedom):
867     case(Nodes):
868     case(Elements):
869     case(FaceElements):
870     case(Points):
871     case(ContactElementsZero):
872     case(ContactElementsOne):
873     return true;
874     default:
875     Finley_ErrorCode=TYPE_ERROR;
876     sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);
877     }
878     break;
879     case(ReducedDegreesOfFreedom):
880     switch(functionSpaceType_target) {
881     case(ReducedDegreesOfFreedom):
882     case(Nodes):
883     case(Elements):
884     case(FaceElements):
885     case(Points):
886     case(ContactElementsZero):
887     case(ContactElementsOne):
888     return true;
889     case(DegreesOfFreedom):
890     return false;
891     default:
892     Finley_ErrorCode=TYPE_ERROR;
893     sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_target);
894     }
895     break;
896     default:
897     Finley_ErrorCode=TYPE_ERROR;
898     sprintf(Finley_ErrorMsg,"Interpolation On Domain: Finley does not know anything about function space type %d",functionSpaceType_source);
899     break;
900     }
901     checkFinleyError();
902     return false;
903     }
904     bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
905     {
906     return false;
907     }
908     bool MeshAdapter::operator==(const MeshAdapter& other) const
909     {
910     return (m_finleyMesh==other.m_finleyMesh);
911     }
912    
913 jgs 102 int MeshAdapter::getSystemMatrixTypeId(const int solver, const bool symmetry) const
914     {
915     int out=Finley_SystemMatrix_getSystemMatrixTypeId(solver,symmetry?1:0);
916     checkFinleyError();
917     return out;
918     }
919     Data MeshAdapter::getX() const
920     {
921     return continuousFunction(asAbstractContinuousDomain()).getX();
922     }
923     Data MeshAdapter::getNormal() const
924     {
925     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
926     }
927     Data MeshAdapter::getSize() const
928     {
929     return function(asAbstractContinuousDomain()).getSize();
930     }
931    
932    
933 jgs 82 bool MeshAdapter::operator!=(const MeshAdapter& other) const
934     {
935     return !operator==(other);
936     }
937     // bool MeshAdapter::operator==(const AbstractDomain& other) const
938     // {
939     // try {
940     // const MeshAdapter& ref = dynamic_cast<const MeshAdapter&>(other);
941     // return (operator==(ref));
942     // }
943     // catch (bad_cast) {
944     // return false;
945     // }
946     // }
947     } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26