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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (13 years, 7 months ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 38207 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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