/[escript]/trunk/ripley/src/RipleyDomain.cpp
ViewVC logotype

Annotation of /trunk/ripley/src/RipleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3670 - (hide annotations)
Wed Nov 16 02:21:18 2011 UTC (7 years, 10 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp
File size: 80879 byte(s)
First steps...

1 caltinay 3670
2     /*******************************************************
3     *
4     * Copyright (c) 2003-2011 by University of Queensland
5     * Earth Systems Science Computational Center (ESSCC)
6     * http://www.uq.edu.au/esscc
7     *
8     * Primary Business: Queensland, Australia
9     * Licensed under the Open Software License version 3.0
10     * http://www.opensource.org/licenses/osl-3.0.php
11     *
12     *******************************************************/
13    
14     #include <ripley/RipleyDomain.h>
15     #include <ripley/IndexList.h>
16     #include <ripley/Util.h>
17     extern "C" {
18     #include "esysUtils/blocktimer.h"
19     }
20     #include <escript/Data.h>
21     #include <escript/DataFactory.h>
22    
23     #ifdef USE_NETCDF
24     #include <netcdfcpp.h>
25     #endif
26    
27     #ifdef USE_PARMETIS
28     #include <parmetis.h>
29     #endif
30    
31     #include <boost/python/import.hpp>
32     #include <boost/python/tuple.hpp>
33     #include <iomanip>
34    
35     using namespace std;
36    
37     namespace ripley {
38    
39     //
40     // define the static constants
41     RipleyDomain::FunctionSpaceNamesMapType RipleyDomain::m_functionSpaceTypeNames;
42    
43     RipleyDomain::RipleyDomain(const string& name, dim_t numDim, Esys_MPIInfo* mpiInfo) :
44     m_name(name)
45     {
46     setFunctionSpaceTypeNames();
47     m_mpiInfo = Esys_MPIInfo_getReference(mpiInfo);
48     m_nodes.reset(new NodeFile(numDim, m_mpiInfo));
49     if (numDim==3) {
50     m_elements.reset(new ElementFile(Hex8, m_mpiInfo));
51     m_faceElements.reset(new ElementFile(Rec4, m_mpiInfo));
52     } else if (numDim==2) {
53     m_elements.reset(new ElementFile(Rec4, m_mpiInfo));
54     m_faceElements.reset(new ElementFile(Line2, m_mpiInfo));
55     } else
56     throw RipleyException("RipleyDomain: Invalid number of dimensions");
57     m_points.reset(new ElementFile(Point1, m_mpiInfo));
58     m_fullFullPattern = NULL;
59     m_fullReducedPattern = NULL;
60     m_reducedFullPattern = NULL;
61     m_reducedReducedPattern = NULL;
62     }
63    
64     RipleyDomain::~RipleyDomain()
65     {
66     Esys_MPIInfo_free(m_mpiInfo);
67     }
68    
69     int RipleyDomain::getMPISize() const
70     {
71     return m_mpiInfo->size;
72     }
73    
74     int RipleyDomain::getMPIRank() const
75     {
76     return m_mpiInfo->rank;
77     }
78    
79     void RipleyDomain::MPIBarrier() const
80     {
81     #ifdef ESYS_MPI
82     MPI_Barrier(m_mpiInfo->comm);
83     #endif
84     }
85    
86     bool RipleyDomain::onMasterProcessor() const
87     {
88     return m_mpiInfo->rank == 0;
89     }
90    
91     #ifdef ESYS_MPI
92     MPI_Comm
93     #else
94     unsigned int
95     #endif
96     RipleyDomain::getMPIComm() const
97     {
98     #ifdef ESYS_MPI
99     return m_mpiInfo->comm;
100     #else
101     return 0;
102     #endif
103     }
104    
105     void RipleyDomain::write(const string& filename) const
106     {
107     if (m_mpiInfo->size > 1)
108     throw RipleyException("write: only single processor runs are supported.");
109    
110     // open file
111     ofstream f(filename.c_str());
112     if (!f.is_open()) {
113     stringstream msg;
114     msg << "write: Opening file " << filename << " for writing failed.";
115     throw RipleyException(msg.str());
116     }
117    
118     // write header
119     f << m_name << endl;
120    
121     f.precision(15);
122     f.setf(ios::scientific, ios::floatfield);
123    
124     // write nodes
125     dim_t numDim = getDim();
126     f << numDim << "D-Nodes " << m_nodes->getNumNodes() << endl;
127     for (dim_t i = 0; i < m_nodes->getNumNodes(); i++) {
128     f << m_nodes->getIdVector()[i] << " "
129     << m_nodes->getGlobalDegreesOfFreedom()[i] << " "
130     << m_nodes->getTagVector()[i];
131     for (dim_t j = 0; j < numDim; j++)
132     f << " " << setw(20) << m_nodes->getCoordinates()[INDEX2(j, i, numDim)];
133     f << endl;
134     }
135    
136     // write all element types
137     const ElementFile_ptr eFiles[3] = { m_elements, m_faceElements, m_points };
138    
139     for (size_t k=0; k<3; k++) {
140     f << eFiles[k]->getName() << " " << eFiles[k]->getNumElements() << endl;
141     dim_t NN = eFiles[k]->getNumNodes();
142     for (dim_t i = 0; i < eFiles[k]->getNumElements(); i++) {
143     f << eFiles[k]->getIdVector()[i] << " "
144     << eFiles[k]->getTagVector()[i];
145     for (dim_t j = 0; j < NN; j++)
146     f << " " << m_nodes->getIdVector()[eFiles[k]->getNodes()[INDEX2(j, i, NN)]];
147     f << endl;
148     }
149     }
150    
151     // write tags
152     if (m_tagMap.size() > 0) {
153     f << "Tags" << endl;
154     TagMap::const_iterator it;
155     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
156     f << it->first << " " << it->second << endl;
157     }
158     }
159     f.close();
160     #ifdef Ripley_TRACE
161     cout << "mesh " << m_name << " has been written to file " << filename << endl;
162     #endif
163     }
164    
165     void RipleyDomain::Print_Mesh_Info(const bool full) const
166     {
167     cout << "Ripley_PrintMesh_Info running on CPU " <<
168     m_mpiInfo->rank << " of " << m_mpiInfo->size << endl;
169     cout << "\tMesh name '" << m_name << "'" << endl;
170    
171     // write nodes
172     int numDim = getDim();
173     cout << "\tNodes: " << numDim << "D-Nodes " << m_nodes->getNumNodes() << endl;
174     if (full) {
175     cout << "\t Id Tag gDOF gNI grDfI grNI: Coordinates" << endl;
176     cout.precision(15);
177     cout.setf(ios::scientific, ios::floatfield);
178     for (int i = 0; i < m_nodes->getNumNodes(); i++) {
179     cout << "\t " << setw(5) << m_nodes->getIdVector()[i] << " "
180     << setw(5) << m_nodes->getTagVector()[i] << " "
181     << setw(5) << m_nodes->getGlobalDegreesOfFreedom()[i] << " "
182     << setw(5) << m_nodes->getGlobalNodesIndex()[i] << " "
183     << setw(5) << m_nodes->getGlobalReducedDOFIndex()[i] << " "
184     << setw(5) << m_nodes->getGlobalReducedNodesIndex()[i] << ": ";
185     for (int j = 0; j < numDim; j++)
186     cout << " " << m_nodes->getCoordinates()[INDEX2(j, i, numDim)];
187     cout << endl;
188     }
189     }
190    
191     // write all element types
192     const char *eNames[3] = { "Elements", "Face elements", "Points" };
193     const ElementFile_ptr eFiles[3] = { m_elements, m_faceElements, m_points };
194    
195     for (size_t k=0; k<3; k++) {
196     int mine = 0, overlap = 0;
197     for (dim_t i = 0; i < eFiles[k]->getNumElements(); i++) {
198     if (eFiles[k]->getOwnerVector()[i] == m_mpiInfo->rank)
199     mine++;
200     else
201     overlap++;
202     }
203     cout << "\t" << eNames[k] << ": " << eFiles[k]->getName()
204     << " " << eFiles[k]->getNumElements()
205     << " (TypeId=" << eFiles[k]->getTypeId() << ") owner="
206     << mine << " overlap=" << overlap << endl;
207     int NN = eFiles[k]->getNumNodes();
208     if (full) {
209     cout << "\t Id Tag Owner Color: Nodes" << endl;
210     for (int i = 0; i < eFiles[k]->getNumElements(); i++) {
211     cout << "\t " << setw(5) << eFiles[k]->getIdVector()[i] << " "
212     << setw(5) << eFiles[k]->getTagVector()[i] << " "
213     << setw(5) << eFiles[k]->getOwnerVector()[i] << " "
214     << setw(5) << eFiles[k]->getColorVector()[i] << ": ";
215     for (int j = 0; j < NN; j++)
216     cout << " " << setw(5) << m_nodes->getIdVector()[eFiles[k]->getNodes()[INDEX2(j, i, NN)]];
217     cout << endl;
218     }
219     }
220     }
221    
222    
223     // write tags
224     if (m_tagMap.size() > 0) {
225     cout << "\tTags:" << endl;
226     TagMap::const_iterator it;
227     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
228     cout << "\t " << setw(5) << it->second << " "
229     << it->first << endl;
230     }
231     }
232     }
233    
234     void RipleyDomain::dump(const string& fileName) const
235     {
236     #ifdef USE_NETCDF
237     NcDim* ncdim = NULL;
238     int mpi_size = m_mpiInfo->size;
239     int mpi_rank = m_mpiInfo->rank;
240     int numDim = getDim();
241    
242     /* Incoming token indicates it's my turn to write */
243     #ifdef ESYS_MPI
244     MPI_Status status;
245     if (mpi_rank>0) {
246     int dummy;
247     MPI_Recv(&dummy, 0, MPI_INT, mpi_rank-1, 81800, m_mpiInfo->comm, &status);
248     }
249     #endif
250    
251     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
252     mpi_size, mpi_rank);
253    
254     // NetCDF error handler
255     NcError err(NcError::verbose_nonfatal);
256     // Create the file.
257     NcFile dataFile(newFileName, NcFile::Replace);
258     const string msgPrefix("dump: netCDF operation failed - ");
259    
260     // check if writing was successful
261     if (!dataFile.is_valid())
262     throw RipleyException(msgPrefix+"Open file for output");
263    
264     const size_t numTags = m_tagMap.size();
265    
266     if (numTags > 0)
267     if (! (ncdim = dataFile.add_dim("dim_Tags", numTags)) )
268     throw RipleyException(msgPrefix+"add_dim(dim_Tags)");
269    
270     // Attributes: MPI size, MPI rank, Name, order, reduced_order
271     if (!dataFile.add_att("mpi_size", mpi_size))
272     throw RipleyException(msgPrefix+"add_att(mpi_size)");
273     if (!dataFile.add_att("mpi_rank", mpi_rank))
274     throw RipleyException(msgPrefix+"add_att(mpi_rank)");
275     if (!dataFile.add_att("Name", m_name.c_str()))
276     throw RipleyException(msgPrefix+"add_att(Name)");
277     if (!dataFile.add_att("numDim", numDim))
278     throw RipleyException(msgPrefix+"add_att(numDim)");
279     if (!dataFile.add_att("order", 1))
280     throw RipleyException(msgPrefix+"add_att(order)");
281     if (!dataFile.add_att("reduced_order", 1))
282     throw RipleyException(msgPrefix+"add_att(reduced_order)");
283     if (!dataFile.add_att("num_Tags", static_cast<int>(numTags)))
284     throw RipleyException(msgPrefix+"add_att(num_Tags)");
285    
286     m_nodes->dumpToNetCDF(dataFile);
287     m_elements->dumpToNetCDF(dataFile, "Elements");
288     m_faceElements->dumpToNetCDF(dataFile, "FaceElements");
289     m_points->dumpToNetCDF(dataFile, "Points");
290    
291     // // // // // TagMap // // // // //
292     if (numTags > 0) {
293    
294     // Copy tag keys into temp array
295     vector<int> Tags_keys;
296     TagMap::const_iterator it;
297     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
298     Tags_keys.push_back(it->second);
299     }
300    
301     // Tags_keys
302     NcVar *ncVar;
303     if (! (ncVar = dataFile.add_var("Tags_keys", ncInt, ncdim)) )
304     throw RipleyException(msgPrefix+"add_var(Tags_keys)");
305     if (!ncVar->put(&Tags_keys[0], numTags))
306     throw RipleyException(msgPrefix+"put(Tags_keys)");
307    
308     // Tags_names_*
309     // This is an array of strings, it should be stored as an array but
310     // instead I have hacked in one attribute per string because the NetCDF
311     // manual doesn't tell how to do an array of strings
312     int i = 0;
313     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
314     stringstream name;
315     name << "Tags_name_" << i;
316     if (!dataFile.add_att(name.str().c_str(), it->first.c_str()) )
317     throw RipleyException(msgPrefix+"add_att(Tags_names_XX)");
318     i++;
319     }
320     }
321    
322     // Send token to next MPI process so he can take his turn
323     #ifdef ESYS_MPI
324     if (mpi_rank<mpi_size-1) {
325     int dummy = 0;
326     MPI_Send(&dummy, 0, MPI_INT, mpi_rank+1, 81800, m_mpiInfo->comm);
327     }
328     #endif
329    
330     // netCDF file is closed by destructor of NcFile object
331     #else
332     throw RipleyException("dump: not configured with netCDF. Please contact your installation manager.");
333     #endif /* USE_NETCDF */
334     }
335    
336     string RipleyDomain::getDescription() const
337     {
338     return "RipleyMesh";
339     }
340    
341     string RipleyDomain::functionSpaceTypeAsString(int functionSpaceType) const
342     {
343     FunctionSpaceNamesMapType::iterator loc;
344     loc=m_functionSpaceTypeNames.find(functionSpaceType);
345     if (loc!=m_functionSpaceTypeNames.end())
346     return loc->second;
347    
348     return "Invalid function space type code";
349     }
350    
351     bool RipleyDomain::isValidFunctionSpaceType(int functionSpaceType) const
352     {
353     FunctionSpaceNamesMapType::iterator loc;
354     loc=m_functionSpaceTypeNames.find(functionSpaceType);
355     return (loc!=m_functionSpaceTypeNames.end());
356     }
357    
358     void RipleyDomain::setFunctionSpaceTypeNames()
359     {
360     if (m_functionSpaceTypeNames.size() > 0)
361     return;
362    
363     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Ripley_DegreesOfFreedom"));
364     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Ripley_ReducedDegreesOfFreedom"));
365     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(Nodes,"Ripley_Nodes"));
366     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedNodes,"Ripley_Reduced_Nodes"));
367     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(Elements,"Ripley_Elements"));
368     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedElements,"Ripley_Reduced_Elements"));
369     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(FaceElements,"Ripley_Face_Elements"));
370     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Ripley_Reduced_Face_Elements"));
371     m_functionSpaceTypeNames.insert(FunctionSpaceNamesMapType::value_type(Points,"Ripley_Points"));
372     }
373    
374     int RipleyDomain::getContinuousFunctionCode() const
375     {
376     return Nodes;
377     }
378    
379     int RipleyDomain::getReducedContinuousFunctionCode() const
380     {
381     return ReducedNodes;
382     }
383    
384     int RipleyDomain::getFunctionCode() const
385     {
386     return Elements;
387     }
388    
389     int RipleyDomain::getReducedFunctionCode() const
390     {
391     return ReducedElements;
392     }
393    
394     int RipleyDomain::getFunctionOnBoundaryCode() const
395     {
396     return FaceElements;
397     }
398    
399     int RipleyDomain::getReducedFunctionOnBoundaryCode() const
400     {
401     return ReducedFaceElements;
402     }
403    
404     int RipleyDomain::getFunctionOnContactZeroCode() const
405     {
406     throw RipleyException("Ripley does not support contact elements");
407     }
408    
409     int RipleyDomain::getReducedFunctionOnContactZeroCode() const
410     {
411     throw RipleyException("Ripley does not support contact elements");
412     }
413    
414     int RipleyDomain::getFunctionOnContactOneCode() const
415     {
416     throw RipleyException("Ripley does not support contact elements");
417     }
418    
419     int RipleyDomain::getReducedFunctionOnContactOneCode() const
420     {
421     throw RipleyException("Ripley does not support contact elements");
422     }
423    
424     int RipleyDomain::getSolutionCode() const
425     {
426     return DegreesOfFreedom;
427     }
428    
429     int RipleyDomain::getReducedSolutionCode() const
430     {
431     return ReducedDegreesOfFreedom;
432     }
433    
434     int RipleyDomain::getDiracDeltaFunctionsCode() const
435     {
436     return Points;
437     }
438    
439     //
440     // returns the spatial dimension of the mesh
441     //
442     int RipleyDomain::getDim() const
443     {
444     return m_nodes->getNumDim();
445     }
446    
447     //
448     // returns the number of data points summed across all MPI processes
449     //
450     int RipleyDomain::getNumDataPointsGlobal() const
451     {
452     return m_nodes->getGlobalNumNodes();
453     }
454    
455     //
456     // returns the number of data points per sample and the number of samples
457     // needed to represent data on a parts of the mesh.
458     //
459     pair<int,int> RipleyDomain::getDataShape(int functionSpaceCode) const
460     {
461     int numDataPointsPerSample=0;
462     int numSamples=0;
463     switch (functionSpaceCode) {
464     case Nodes:
465     numDataPointsPerSample=1;
466     numSamples=m_nodes->getNumNodes();
467     break;
468     case ReducedNodes:
469     numDataPointsPerSample=1;
470     numSamples=m_nodes->getNumReducedNodes();
471     break;
472     case Elements:
473     numSamples=m_elements->getNumElements();
474     numDataPointsPerSample=m_elements->getNumLocalDim()+1;
475     break;
476     case ReducedElements:
477     numSamples=m_elements->getNumElements();
478     numDataPointsPerSample=(m_elements->getNumLocalDim()==0)?0:1;
479     break;
480     case FaceElements:
481     numDataPointsPerSample=m_faceElements->getNumLocalDim()+1;
482     numSamples=m_faceElements->getNumElements();
483     break;
484     case ReducedFaceElements:
485     numDataPointsPerSample=(m_faceElements->getNumLocalDim()==0)?0:1;
486     numSamples=m_faceElements->getNumElements();
487     break;
488     case Points:
489     numDataPointsPerSample=1;
490     numSamples=m_points->getNumElements();
491     break;
492     case DegreesOfFreedom:
493     numDataPointsPerSample=1;
494     numSamples=m_nodes->getNumDegreesOfFreedom();
495     break;
496     case ReducedDegreesOfFreedom:
497     numDataPointsPerSample=1;
498     numSamples=m_nodes->getNumReducedDegreesOfFreedom();
499     break;
500     default:
501     stringstream temp;
502     temp << "Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
503     throw RipleyException(temp.str());
504     }
505     return pair<int,int>(numDataPointsPerSample,numSamples);
506     }
507    
508     //
509     // adds linear PDE of second order into a given stiffness matrix and right hand side
510     //
511     void RipleyDomain::addPDEToSystem(
512     escript::AbstractSystemMatrix& mat, escript::Data& rhs,
513     const escript::Data& A, const escript::Data& B, const escript::Data& C,
514     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
515     const escript::Data& d, const escript::Data& y,
516     const escript::Data& d_contact, const escript::Data& y_contact,
517     const escript::Data& d_dirac,const escript::Data& y_dirac) const
518     {
519     if (!d_contact.isEmpty() || !y_contact.isEmpty())
520     throw RipleyException("Ripley does not support contact elements");
521    
522     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
523     if (smat==NULL)
524     throw RipleyException("Ripley only accepts its own system matrices");
525     /*
526     escriptDataC _rhs=rhs.getDataC();
527     escriptDataC _A =A.getDataC();
528     escriptDataC _B=B.getDataC();
529     escriptDataC _C=C.getDataC();
530     escriptDataC _D=D.getDataC();
531     escriptDataC _X=X.getDataC();
532     escriptDataC _Y=Y.getDataC();
533     escriptDataC _d=d.getDataC();
534     escriptDataC _y=y.getDataC();
535     escriptDataC _d_dirac=d_dirac.getDataC();
536     escriptDataC _y_dirac=y_dirac.getDataC();
537    
538     Ripley_Mesh* mesh=m_ripleyMesh.get();
539     Ripley_Assemble_PDE(m_nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y);
540     checkPasoError();
541    
542     Ripley_Assemble_PDE(m_nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y);
543     checkPasoError();
544    
545     Ripley_Assemble_PDE(m_nodes,mesh->Points, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d_dirac, 0, &_y_dirac);
546     checkPasoError();
547     */
548     }
549    
550     void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
551     const escript::Data& D, const escript::Data& d,
552     const escript::Data& d_dirac, const bool useHRZ) const
553     {
554     /*
555     escriptDataC _mat=mat.getDataC();
556     escriptDataC _D=D.getDataC();
557     escriptDataC _d=d.getDataC();
558     escriptDataC _d_dirac=d_dirac.getDataC();
559    
560     Ripley_Mesh* mesh=m_ripleyMesh.get();
561     Ripley_Assemble_LumpedSystem(m_nodes, mesh->Elements, &_mat, &_D, useHRZ);
562     checkPasoError();
563    
564     Ripley_Assemble_LumpedSystem(m_nodes, mesh->FaceElements, &_mat, &_d, useHRZ);
565     checkPasoError();
566    
567     Ripley_Assemble_LumpedSystem(m_nodes, mesh->FaceElements, &_mat, &_d_dirac, useHRZ);
568     checkPasoError();
569     */
570     }
571    
572    
573     //
574     // adds linear PDE of second order into the right hand side only
575     //
576     void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
577     const escript::Data& Y, const escript::Data& y,
578     const escript::Data& y_contact, const escript::Data& y_dirac) const
579     {
580     if (!y_contact.isEmpty())
581     throw RipleyException("Ripley does not support contact elements");
582     /*
583     Ripley_Mesh* mesh=m_ripleyMesh.get();
584    
585     escriptDataC _rhs=rhs.getDataC();
586     escriptDataC _X=X.getDataC();
587     escriptDataC _Y=Y.getDataC();
588     escriptDataC _y=y.getDataC();
589     escriptDataC _y_dirac=y_dirac.getDataC();
590     Ripley_Assemble_PDE(m_nodes, mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
591     checkPasoError();
592    
593     Ripley_Assemble_PDE(m_nodes, mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
594     checkPasoError();
595    
596     Ripley_Assemble_PDE(m_nodes, mesh->Points, 0, &_rhs, 0, 0, 0, 0, 0, &_y_dirac );
597     checkPasoError();
598     */
599     }
600    
601     //
602     // adds PDE of second order into a transport problem
603     //
604     void RipleyDomain::addPDEToTransportProblem(
605     escript::AbstractTransportProblem& tp,
606     escript::Data& source, const escript::Data& M,
607     const escript::Data& A, const escript::Data& B, const escript::Data& C,
608     const escript::Data& D, const escript::Data& X, const escript::Data& Y,
609     const escript::Data& d, const escript::Data& y,
610     const escript::Data& d_contact, const escript::Data& y_contact,
611     const escript::Data& d_dirac, const escript::Data& y_dirac) const
612     {
613     if (!d_contact.isEmpty() || !y_contact.isEmpty())
614     throw RipleyException("Ripley does not support contact elements");
615     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
616     if (tpa==NULL)
617     throw RipleyException("Ripley only accepts its own transport problems");
618     /*
619     DataTypes::ShapeType shape;
620     source.expand();
621     escriptDataC _source=source.getDataC();
622     escriptDataC _M=M.getDataC();
623     escriptDataC _A=A.getDataC();
624     escriptDataC _B=B.getDataC();
625     escriptDataC _C=C.getDataC();
626     escriptDataC _D=D.getDataC();
627     escriptDataC _X=X.getDataC();
628     escriptDataC _Y=Y.getDataC();
629     escriptDataC _d=d.getDataC();
630     escriptDataC _y=y.getDataC();
631     escriptDataC _d_dirac=d_dirac.getDataC();
632     escriptDataC _y_dirac=y_dirac.getDataC();
633    
634     Ripley_Mesh* mesh=m_ripleyMesh.get();
635     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
636     Ripley_Assemble_PDE(m_nodes, mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0);
637     checkPasoError();
638    
639     Ripley_Assemble_PDE(m_nodes, mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y);
640     checkPasoError();
641    
642     Ripley_Assemble_PDE(m_nodes, mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y);
643     checkPasoError();
644    
645     Ripley_Assemble_PDE(m_nodes, mesh->Points, _tp->transport_matrix, &_source, 0, 0, 0, &_d_dirac, 0, &_y_dirac);
646     checkPasoError();
647     */
648     }
649    
650     //
651     // interpolates data between different function spaces
652     //
653     void RipleyDomain::interpolateOnDomain(escript::Data& target,
654     const escript::Data& in) const
655     {
656     const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
657     const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
658     if (inDomain != *this)
659     throw RipleyException("Illegal domain of interpolant");
660     if (targetDomain != *this)
661     throw RipleyException("Illegal domain of interpolation target");
662     /*
663     Ripley_Mesh* mesh=m_ripleyMesh.get();
664     escriptDataC _target=target.getDataC();
665     escriptDataC _in=in.getDataC();
666    
667     switch (in.getFunctionSpace().getTypeCode()) {
668     case Nodes:
669     switch (target.getFunctionSpace().getTypeCode()) {
670     case Nodes:
671     case ReducedNodes:
672     case DegreesOfFreedom:
673     case ReducedDegreesOfFreedom:
674     Ripley_Assemble_CopyNodalData(m_nodes, &_target, &_in);
675     break;
676     case Elements:
677     case ReducedElements:
678     Ripley_Assemble_interpolate(m_nodes, mesh->Elements, &_in, &_target);
679     break;
680     case FaceElements:
681     case ReducedFaceElements:
682     Ripley_Assemble_interpolate(m_nodes, mesh->FaceElements, &_in, &_target);
683     break;
684     case Points:
685     Ripley_Assemble_interpolate(m_nodes, mesh->Points, &_in, &_target);
686     break;
687     default:
688     stringstream temp;
689     temp << "Interpolation on Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
690     throw RipleyException(temp.str());
691     }
692     break;
693     case ReducedNodes:
694     switch (target.getFunctionSpace().getTypeCode()) {
695     case Nodes:
696     case ReducedNodes:
697     case DegreesOfFreedom:
698     case ReducedDegreesOfFreedom:
699     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);
700     break;
701     case Elements:
702     case ReducedElements:
703     Ripley_Assemble_interpolate(m_nodes,mesh->Elements,&_in,&_target);
704     break;
705     case FaceElements:
706     case ReducedFaceElements:
707     Ripley_Assemble_interpolate(m_nodes,mesh->FaceElements,&_in,&_target);
708     break;
709     case Points:
710     Ripley_Assemble_interpolate(m_nodes,mesh->Points,&_in,&_target);
711     break;
712     default:
713     stringstream temp;
714     temp << "Interpolation on Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
715     throw RipleyException(temp.str());
716     }
717     break;
718     case Elements:
719     if (target.getFunctionSpace().getTypeCode()==Elements) {
720     Ripley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
721     } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
722     Ripley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
723     } else {
724     throw RipleyException("No interpolation with data on elements possible.");
725     }
726     break;
727     case ReducedElements:
728     if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
729     Ripley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
730     } else {
731     throw RipleyException("No interpolation with data on elements with reduced integration order possible.");
732     }
733     break;
734     case FaceElements:
735     if (target.getFunctionSpace().getTypeCode()==FaceElements) {
736     Ripley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
737     } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
738     Ripley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
739     } else {
740     throw RipleyException("No interpolation with data on face elements possible.");
741     }
742     break;
743     case ReducedFaceElements:
744     if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
745     Ripley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
746     } else {
747     throw RipleyException("No interpolation with data on face elements with reduced integration order possible.");
748     }
749     break;
750     case Points:
751     if (target.getFunctionSpace().getTypeCode()==Points) {
752     Ripley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
753     } else {
754     throw RipleyException("No interpolation with data on points possible.");
755     }
756     break;
757     case DegreesOfFreedom:
758     switch (target.getFunctionSpace().getTypeCode()) {
759     case ReducedDegreesOfFreedom:
760     case DegreesOfFreedom:
761     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);
762     break;
763    
764     case Nodes:
765     case ReducedNodes:
766     if (getMPISize()>1) {
767     escript::Data temp=escript::Data(in);
768     temp.expand();
769     escriptDataC _in2 = temp.getDataC();
770     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in2);
771     } else {
772     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);
773     }
774     break;
775     case Elements:
776     case ReducedElements:
777     if (getMPISize()>1) {
778     escript::Data temp=escript::Data(in, continuousFunction(*this) );
779     escriptDataC _in2 = temp.getDataC();
780     Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in2,&_target);
781     } else {
782     Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in,&_target);
783     }
784     break;
785     case FaceElements:
786     case ReducedFaceElements:
787     if (getMPISize()>1) {
788     escript::Data temp=escript::Data(in, continuousFunction(*this) );
789     escriptDataC _in2 = temp.getDataC();
790     Ripley_Assemble_interpolate(m_nodes, mesh->FaceElements,&_in2,&_target);
791     } else {
792     Ripley_Assemble_interpolate(m_nodes, mesh->FaceElements,&_in,&_target);
793     }
794     break;
795     case Points:
796     if (getMPISize()>1) {
797     //escript::Data temp=escript::Data(in, continuousFunction(*this) );
798     //escriptDataC _in2 = temp.getDataC();
799     } else {
800     Ripley_Assemble_interpolate(m_nodes,mesh->Points,&_in,&_target);
801     }
802     break;
803     default:
804     stringstream temp;
805     temp << "Interpolation On Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
806     throw RipleyException(temp.str());
807     }
808     break;
809     case ReducedDegreesOfFreedom:
810     switch (target.getFunctionSpace().getTypeCode()) {
811     case Nodes:
812     throw RipleyException("Ripley does not support interpolation from reduced degrees of freedom to mesh nodes.");
813     break;
814     case ReducedNodes:
815     if (getMPISize()>1) {
816     escript::Data temp=escript::Data(in);
817     temp.expand();
818     escriptDataC _in2 = temp.getDataC();
819     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in2);
820     } else {
821     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);
822     }
823     break;
824     case DegreesOfFreedom:
825     throw RipleyException("Ripley does not support interpolation from reduced degrees of freedom to degrees of freedom");
826     case ReducedDegreesOfFreedom:
827     Ripley_Assemble_CopyNodalData(m_nodes,&_target,&_in);
828     break;
829     case Elements:
830     case ReducedElements:
831     if (getMPISize()>1) {
832     escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );
833     escriptDataC _in2 = temp.getDataC();
834     Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in2,&_target);
835     } else {
836     Ripley_Assemble_interpolate(m_nodes, mesh->Elements,&_in,&_target);
837     }
838     break;
839     case FaceElements:
840     case ReducedFaceElements:
841     if (getMPISize()>1) {
842     escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );
843     escriptDataC _in2 = temp.getDataC();
844     Ripley_Assemble_interpolate(m_nodes,mesh->FaceElements,&_in2,&_target);
845     } else {
846     Ripley_Assemble_interpolate(m_nodes,mesh->FaceElements,&_in,&_target);
847     }
848     break;
849     case Points:
850     if (getMPISize()>1) {
851     escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );
852     escriptDataC _in2 = temp.getDataC();
853     Ripley_Assemble_interpolate(m_nodes, mesh->Points,&_in2,&_target);
854     } else {
855     Ripley_Assemble_interpolate(m_nodes, mesh->Points,&_in,&_target);
856     }
857     break;
858     default:
859     stringstream temp;
860     temp << "Interpolation On Domain: Ripley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
861     throw RipleyException(temp.str());
862     }
863     break;
864     default:
865     stringstream temp;
866     temp << "Interpolation On Domain: Ripley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
867     throw RipleyException(temp.str());
868     }
869     */
870     checkPasoError();
871     }
872    
873     //
874     // copies the locations of sample points into arg:
875     //
876     void RipleyDomain::setToX(escript::Data& arg) const
877     {
878     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(*(arg.getFunctionSpace().getDomain()));
879     if (argDomain != *this)
880     throw RipleyException("setToX: Illegal domain of data point locations");
881    
882     // do we need to interpolate?
883     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
884     m_nodes->assembleCoordinates(arg);
885     } else {
886     escript::Data tmp_data=Vector(0.0, continuousFunction(*this), true);
887     m_nodes->assembleCoordinates(tmp_data);
888     // this is then interpolated onto arg:
889     interpolateOnDomain(arg, tmp_data);
890     }
891     }
892    
893     //
894     // returns the normal vectors at the location of data points as a Data object
895     //
896     void RipleyDomain::setToNormal(escript::Data& normal) const
897     {
898     /*
899     const RipleyDomain& normalDomain=dynamic_cast<const RipleyDomain&>(*(normal.getFunctionSpace().getDomain()));
900     if (normalDomain!=*this)
901     throw RipleyException("Illegal domain of normal locations");
902     Ripley_Mesh* mesh=m_ripleyMesh.get();
903     escriptDataC _normal=normal.getDataC();
904     switch(normal.getFunctionSpace().getTypeCode()) {
905     case Nodes:
906     throw RipleyException("Ripley does not support surface normal vectors for nodes");
907     case ReducedNodes:
908     throw RipleyException("Ripley does not support surface normal vectors for reduced nodes");
909     case Elements:
910     throw RipleyException("Ripley does not support surface normal vectors for elements");
911     case ReducedElements:
912     throw RipleyException("Ripley does not support surface normal vectors for elements with reduced integration order");
913     case FaceElements:
914     Ripley_Assemble_setNormal(m_nodes,mesh->FaceElements,&_normal);
915     break;
916     case ReducedFaceElements:
917     Ripley_Assemble_setNormal(m_nodes,mesh->FaceElements,&_normal);
918     break;
919     case Points:
920     throw RipleyException("Ripley does not support surface normal vectors for point elements");
921     case DegreesOfFreedom:
922     throw RipleyException("Ripley does not support surface normal vectors for degrees of freedom.");
923     case ReducedDegreesOfFreedom:
924     throw RipleyException("Ripley does not support surface normal vectors for reduced degrees of freedom.");
925     default:
926     stringstream temp;
927     temp << "Normal Vectors: Ripley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
928     throw RipleyException(temp.str());
929     }
930     */
931     }
932    
933     //
934     // interpolates data to other domain
935     //
936     void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
937     {
938     escript::const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
939     const RipleyDomain* targetDomain=dynamic_cast<const RipleyDomain*>(targetDomain_p.get());
940     if (targetDomain!=this)
941     throw RipleyException("Illegal domain of interpolation target");
942    
943     throw RipleyException("Ripley does not allow interpolation across domains yet");
944     }
945    
946     //
947     // calculates the integral of a function arg
948     //
949     void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
950     {
951     /*
952     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(*(arg.getFunctionSpace().getDomain()));
953     if (argDomain!=*this)
954     throw RipleyException("Illegal domain of integration kernel");
955    
956     double blocktimer_start = blocktimer_time();
957     Ripley_Mesh* mesh=m_ripleyMesh.get();
958     escriptDataC _temp;
959     escript::Data temp;
960     escriptDataC _arg=arg.getDataC();
961     switch(arg.getFunctionSpace().getTypeCode()) {
962     case Nodes:
963     case ReducedNodes:
964     case DegreesOfFreedom:
965     case ReducedDegreesOfFreedom:
966     temp=escript::Data( arg, escript::function(*this) );
967     _temp=temp.getDataC();
968     Ripley_Assemble_integrate(m_nodes,mesh->Elements,&_temp,&integrals[0]);
969     break;
970     case Elements:
971     case ReducedElements:
972     Ripley_Assemble_integrate(m_nodes,mesh->Elements,&_arg,&integrals[0]);
973     break;
974     case FaceElements:
975     case ReducedFaceElements:
976     Ripley_Assemble_integrate(m_nodes,mesh->FaceElements,&_arg,&integrals[0]);
977     break;
978     case Points:
979     throw RipleyException("Integral of data on points is not supported.");
980     default:
981     stringstream temp;
982     temp << "Integrals: Ripley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
983     throw RipleyException(temp.str());
984     }
985     blocktimer_increment("integrate()", blocktimer_start);
986     */
987     }
988    
989     //
990     // calculates the gradient of arg
991     //
992     void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
993     {
994     /*
995     const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(*(arg.getFunctionSpace().getDomain()));
996     if (argDomain!=*this)
997     throw RipleyException("Illegal domain of gradient argument");
998     const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(*(grad.getFunctionSpace().getDomain()));
999     if (gradDomain!=*this)
1000     throw RipleyException("Illegal domain of gradient");
1001    
1002     Ripley_Mesh* mesh=m_ripleyMesh.get();
1003     escriptDataC _grad=grad.getDataC();
1004     escriptDataC nodeDataC;
1005     escript::Data temp;
1006     if (getMPISize()>1) {
1007     if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1008     temp=escript::Data(arg, continuousFunction(*this) );
1009     nodeDataC = temp.getDataC();
1010     } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1011     temp=escript::Data(arg, reducedContinuousFunction(*this) );
1012     nodeDataC = temp.getDataC();
1013     } else {
1014     nodeDataC = arg.getDataC();
1015     }
1016     } else {
1017     nodeDataC = arg.getDataC();
1018     }
1019     switch(grad.getFunctionSpace().getTypeCode()) {
1020     case Nodes:
1021     throw RipleyException("Gradient at nodes is not supported.");
1022     case ReducedNodes:
1023     throw RipleyException("Gradient at reduced nodes is not supported.");
1024     case Elements:
1025     Ripley_Assemble_gradient(m_nodes,mesh->Elements,&_grad,&nodeDataC);
1026     break;
1027     case ReducedElements:
1028     Ripley_Assemble_gradient(m_nodes,mesh->Elements,&_grad,&nodeDataC);
1029     break;
1030     case FaceElements:
1031     Ripley_Assemble_gradient(m_nodes,mesh->FaceElements,&_grad,&nodeDataC);
1032     break;
1033     case ReducedFaceElements:
1034     Ripley_Assemble_gradient(m_nodes,mesh->FaceElements,&_grad,&nodeDataC);
1035     break;
1036     case Points:
1037     throw RipleyException("Gradient at points is not supported");
1038     case DegreesOfFreedom:
1039     throw RipleyException("Gradient at degrees of freedom is not supported");
1040     case ReducedDegreesOfFreedom:
1041     throw RipleyException("Gradient at reduced degrees of freedom is not supported");
1042     default:
1043     stringstream temp;
1044     temp << "Gradient: Ripley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
1045     throw RipleyException(temp.str());
1046     }
1047     */
1048     }
1049    
1050     //
1051     // returns the size of elements
1052     //
1053     void RipleyDomain::setToSize(escript::Data& size) const
1054     {
1055     /*
1056     Ripley_Mesh* mesh=m_ripleyMesh.get();
1057     escriptDataC tmp=size.getDataC();
1058     switch(size.getFunctionSpace().getTypeCode()) {
1059     case Nodes:
1060     throw RipleyException("Size of nodes is not supported");
1061     case ReducedNodes:
1062     throw RipleyException("Size of reduced nodes is not supported");
1063     case Elements:
1064     Ripley_Assemble_getSize(m_nodes, mesh->Elements, &tmp);
1065     break;
1066     case ReducedElements:
1067     Ripley_Assemble_getSize(m_nodes, mesh->Elements, &tmp);
1068     break;
1069     case FaceElements:
1070     Ripley_Assemble_getSize(m_nodes, mesh->FaceElements, &tmp);
1071     break;
1072     case ReducedFaceElements:
1073     Ripley_Assemble_getSize(m_nodes, mesh->FaceElements, &tmp);
1074     break;
1075     case Points:
1076     throw RipleyException("Size of point elements is not supported");
1077     case DegreesOfFreedom:
1078     throw RipleyException("Size of degrees of freedom is not supported");
1079     case ReducedDegreesOfFreedom:
1080     throw RipleyException("Size of reduced degrees of freedom is not supported");
1081     default:
1082     stringstream temp;
1083     temp << "Element size: Ripley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
1084     throw RipleyException(temp.str());
1085     }
1086     */
1087     }
1088    
1089     //
1090     // sets the location of nodes
1091     //
1092     void RipleyDomain::setNewX(const escript::Data& new_x)
1093     {
1094     const RipleyDomain& newDomain=dynamic_cast<const RipleyDomain&>(*(new_x.getFunctionSpace().getDomain()));
1095     if (newDomain!=*this)
1096     throw RipleyException("Illegal domain of new node locations");
1097    
1098     if (new_x.getFunctionSpace() == continuousFunction(*this)) {
1099     m_nodes->setCoordinates(new_x);
1100     } else {
1101     escript::Data new_x_inter=escript::Data(new_x, continuousFunction(*this) );
1102     m_nodes->setCoordinates(new_x_inter);
1103     }
1104     }
1105    
1106     bool RipleyDomain::ownSample(int fs_code, index_t id) const
1107     {
1108     #ifdef ESYS_MPI
1109     index_t myFirstNode, myLastNode, k;
1110     if (fs_code == ReducedNodes) {
1111     myFirstNode = m_nodes->getFirstReducedNode();
1112     myLastNode = m_nodes->getLastReducedNode();
1113     k = m_nodes->getIndexForGlobalReducedNode(id);
1114     } else {
1115     myFirstNode = m_nodes->getFirstNode();
1116     myLastNode = m_nodes->getLastNode();
1117     k = m_nodes->getIndexForGlobalNode(id);
1118     }
1119     return static_cast<bool>((myFirstNode <= k) && (k < myLastNode));
1120     #endif
1121     return true;
1122     }
1123    
1124    
1125     //
1126     // creates a SystemMatrixAdapter stiffness matrix and initializes it with zeros
1127     //
1128     escript::ASM_ptr RipleyDomain::newSystemMatrix(const int row_blocksize,
1129     const escript::FunctionSpace& row_functionspace,
1130     const int column_blocksize,
1131     const escript::FunctionSpace& column_functionspace,
1132     const int type) const
1133     {
1134     bool reduceRowOrder=false;
1135     bool reduceColOrder=false;
1136     // is the domain right?
1137     const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
1138     if (row_domain!=*this)
1139     throw RipleyException("Domain of row function space does not match the domain of matrix generator");
1140     const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
1141     if (col_domain!=*this)
1142     throw RipleyException("Domain of columnn function space does not match the domain of matrix generator");
1143     // is the function space type right
1144     if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1145     reduceRowOrder=true;
1146     } else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
1147     throw RipleyException("Illegal function space type for system matrix rows");
1148     if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1149     reduceColOrder=true;
1150     } else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
1151     throw RipleyException("Illegal function space type for system matrix columns");
1152    
1153     // generate matrix:
1154     Paso_SystemMatrixPattern* fsystemMatrixPattern=getPattern(reduceRowOrder, reduceColOrder);
1155     checkPasoError();
1156     Paso_SystemMatrix* fsystemMatrix = Paso_SystemMatrix_alloc(type,
1157     fsystemMatrixPattern, row_blocksize, column_blocksize, FALSE);
1158     checkPasoError();
1159     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1160     SystemMatrixAdapter* sma = new SystemMatrixAdapter(fsystemMatrix,
1161     row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1162     return escript::ASM_ptr(sma);
1163     }
1164    
1165     //
1166     // creates a TransportProblemAdapter
1167     //
1168     escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,
1169     const int blocksize, const escript::FunctionSpace& functionspace,
1170     const int type) const
1171     {
1172     int reduceOrder=0;
1173     // is the domain right?
1174     const RipleyDomain& domain=dynamic_cast<const RipleyDomain&>(*(functionspace.getDomain()));
1175     if (domain!=*this)
1176     throw RipleyException("Domain of function space does not match the domain of transport problem generator");
1177     // is the function space type right
1178     if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1179     reduceOrder=1;
1180     } else if (functionspace.getTypeCode()!=DegreesOfFreedom)
1181     throw RipleyException("Illegal function space type for system matrix rows");
1182    
1183     // generate matrix
1184     Paso_SystemMatrixPattern* fsystemMatrixPattern = getPattern(reduceOrder, reduceOrder);
1185     Paso_TransportProblem* transportProblem = Paso_TransportProblem_alloc(
1186     useBackwardEuler, fsystemMatrixPattern, blocksize);
1187     checkPasoError();
1188     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1189     escript::AbstractTransportProblem* atp=new TransportProblemAdapter(
1190     transportProblem, useBackwardEuler, blocksize, functionspace);
1191     return escript::ATP_ptr(atp);
1192     }
1193    
1194     // returns true if data on the function space is considered as being cell
1195     // centered
1196     bool RipleyDomain::isCellOriented(int functionSpaceCode) const
1197     {
1198     switch(functionSpaceCode) {
1199     case Nodes:
1200     case DegreesOfFreedom:
1201     case ReducedDegreesOfFreedom:
1202     return false;
1203     case Elements:
1204     case FaceElements:
1205     case Points:
1206     case ReducedElements:
1207     case ReducedFaceElements:
1208     return true;
1209     default:
1210     stringstream temp;
1211     temp << "Ripley does not know anything about function space type " << functionSpaceCode;
1212     throw RipleyException(temp.str());
1213     }
1214     return false;
1215     }
1216    
1217     bool RipleyDomain::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1218     {
1219     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1220     class 1: DOF <-> Nodes
1221     class 2: ReducedDOF <-> ReducedNodes
1222     class 3: Points
1223     class 4: Elements
1224     class 5: ReducedElements
1225     class 6: FaceElements
1226     class 7: ReducedFaceElements
1227     class 8: ContactElementZero <-> ContactElementOne
1228     class 9: ReducedContactElementZero <-> ReducedContactElementOne
1229    
1230     There is also a set of lines. Interpolation is possible down a line but not between lines.
1231     class 1 and 2 belong to all lines so aren't considered.
1232     line 0: class 3
1233     line 1: class 4,5
1234     line 2: class 6,7
1235     line 3: class 8,9
1236    
1237     For classes with multiple members (eg class 2) we have vars to record if
1238     there is at least one instance. e.g. hasnodes is true if we have at least
1239     one instance of Nodes.
1240     */
1241     if (fs.empty())
1242     return false;
1243     vector<int> hasclass(10);
1244     vector<int> hasline(4);
1245     bool hasnodes=false;
1246     bool hasrednodes=false;
1247     for (int i=0;i<fs.size();++i) {
1248     switch (fs[i]) {
1249     case Nodes: hasnodes=true; // no break is deliberate
1250     case DegreesOfFreedom:
1251     hasclass[1]=1;
1252     break;
1253     case ReducedNodes: hasrednodes=true; // no break is deliberate
1254     case ReducedDegreesOfFreedom:
1255     hasclass[2]=1;
1256     break;
1257     case Points:
1258     hasline[0]=1;
1259     hasclass[3]=1;
1260     break;
1261     case Elements:
1262     hasclass[4]=1;
1263     hasline[1]=1;
1264     break;
1265     case ReducedElements:
1266     hasclass[5]=1;
1267     hasline[1]=1;
1268     break;
1269     case FaceElements:
1270     hasclass[6]=1;
1271     hasline[2]=1;
1272     break;
1273     case ReducedFaceElements:
1274     hasclass[7]=1;
1275     hasline[2]=1;
1276     break;
1277     default:
1278     return false;
1279     }
1280     }
1281     int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1282    
1283     // fail if we have more than one leaf group
1284     // = there are at least two branches we can't interpolate between
1285     if (totlines>1)
1286     return false;
1287     else if (totlines==1) {
1288     // we have points
1289     if (hasline[0]==1) {
1290     resultcode=Points;
1291     } else if (hasline[1]==1) {
1292     if (hasclass[5]==1) {
1293     resultcode=ReducedElements;
1294     } else {
1295     resultcode=Elements;
1296     }
1297     } else if (hasline[2]==1) {
1298     if (hasclass[7]==1) {
1299     resultcode=ReducedFaceElements;
1300     } else {
1301     resultcode=FaceElements;
1302     }
1303     } else { // so we must be in line3
1304     throw RipleyException("Programmer Error - choosing between contact elements - we should never get here");
1305     }
1306     } else { // totlines==0
1307     if (hasclass[2]==1) {
1308     // something from class 2
1309     resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
1310     } else { // something from class 1
1311     resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
1312     }
1313     }
1314     return true;
1315     }
1316    
1317     bool RipleyDomain::probeInterpolationOnDomain(int functionSpaceType_source,
1318     int functionSpaceType_target) const
1319     {
1320     if (functionSpaceType_target != Nodes &&
1321     functionSpaceType_target != ReducedNodes &&
1322     functionSpaceType_target != ReducedDegreesOfFreedom &&
1323     functionSpaceType_target != DegreesOfFreedom &&
1324     functionSpaceType_target != Elements &&
1325     functionSpaceType_target != ReducedElements &&
1326     functionSpaceType_target != FaceElements &&
1327     functionSpaceType_target != ReducedFaceElements &&
1328     functionSpaceType_target != Points) {
1329     stringstream temp;
1330     temp << "Interpolation On Domain: Ripley does not know anything about function space type " << functionSpaceType_target;
1331     throw RipleyException(temp.str());
1332     }
1333    
1334     switch (functionSpaceType_source) {
1335     case Nodes:
1336     case DegreesOfFreedom:
1337     return true;
1338     case ReducedNodes:
1339     case ReducedDegreesOfFreedom:
1340     return (functionSpaceType_target != Nodes &&
1341     functionSpaceType_target != DegreesOfFreedom);
1342     case Elements:
1343     return (functionSpaceType_target==Elements ||
1344     functionSpaceType_target==ReducedElements);
1345     case ReducedElements:
1346     return (functionSpaceType_target==ReducedElements);
1347     case FaceElements:
1348     return (functionSpaceType_target==FaceElements ||
1349     functionSpaceType_target==ReducedFaceElements);
1350     case ReducedFaceElements:
1351     return (functionSpaceType_target==ReducedFaceElements);
1352     case Points:
1353     return (functionSpaceType_target==Points);
1354    
1355     default:
1356     stringstream temp;
1357     temp << "Interpolation On Domain: Ripley does not know anything about function space type " << functionSpaceType_source;
1358     throw RipleyException(temp.str());
1359     }
1360     return false;
1361     }
1362    
1363     bool RipleyDomain::probeInterpolationACross(int functionSpaceType_source,
1364     const AbstractDomain& targetDomain, int functionSpaceType_target) const
1365     {
1366     return false;
1367     }
1368    
1369     bool RipleyDomain::operator==(const AbstractDomain& other) const
1370     {
1371     const RipleyDomain* temp=dynamic_cast<const RipleyDomain*>(&other);
1372     if (temp != NULL) {
1373     return (m_name==temp->m_name && m_nodes==temp->m_nodes &&
1374     m_elements==temp->m_elements &&
1375     m_faceElements==temp->m_faceElements &&
1376     m_points==temp->m_points);
1377     }
1378     return false;
1379     }
1380    
1381     bool RipleyDomain::operator!=(const AbstractDomain& other) const
1382     {
1383     return !(operator==(other));
1384     }
1385    
1386     int RipleyDomain::getSystemMatrixTypeId(const int solver,
1387     const int preconditioner, const int package, const bool symmetry) const
1388     {
1389     int out=Paso_SystemMatrix_getSystemMatrixTypeId(
1390     SystemMatrixAdapter::mapOptionToPaso(solver),
1391     SystemMatrixAdapter::mapOptionToPaso(preconditioner),
1392     SystemMatrixAdapter::mapOptionToPaso(package),
1393     symmetry?1:0, m_mpiInfo);
1394     checkPasoError();
1395     return out;
1396     }
1397    
1398     int RipleyDomain::getTransportTypeId(const int solver, const int preconditioner,
1399     const int package, const bool symmetry) const
1400     {
1401     int out=Paso_TransportProblem_getTypeId(
1402     SystemMatrixAdapter::mapOptionToPaso(solver),
1403     SystemMatrixAdapter::mapOptionToPaso(preconditioner),
1404     SystemMatrixAdapter::mapOptionToPaso(package),
1405     symmetry?1:0, m_mpiInfo);
1406     checkPasoError();
1407     return out;
1408     }
1409    
1410     escript::Data RipleyDomain::getX() const
1411     {
1412     return continuousFunction(*this).getX();
1413     }
1414    
1415     escript::Data RipleyDomain::getNormal() const
1416     {
1417     return functionOnBoundary(*this).getNormal();
1418     }
1419    
1420     escript::Data RipleyDomain::getSize() const
1421     {
1422     return escript::function(*this).getSize();
1423     }
1424    
1425     const int* RipleyDomain::borrowSampleReferenceIDs(int functionSpaceType) const
1426     {
1427     switch (functionSpaceType) {
1428     case Nodes:
1429     return &m_nodes->getIdVector()[0];
1430     case ReducedNodes:
1431     return &m_nodes->getReducedNodesId()[0];
1432     case Elements:
1433     case ReducedElements:
1434     return &m_elements->getIdVector()[0];
1435     case FaceElements:
1436     case ReducedFaceElements:
1437     return &m_faceElements->getIdVector()[0];
1438     case Points:
1439     return &m_points->getIdVector()[0];
1440     case DegreesOfFreedom:
1441     return &m_nodes->getDegreesOfFreedomId()[0];
1442     case ReducedDegreesOfFreedom:
1443     return &m_nodes->getReducedDegreesOfFreedomId()[0];
1444     default:
1445     stringstream msg;
1446     msg << "borrowSampleReferenceIDs: Invalid function space type " << functionSpaceType << " for domain: " << getDescription();
1447     throw RipleyException(msg.str());
1448     }
1449     }
1450    
1451     int RipleyDomain::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1452     {
1453     throw RipleyException("not implemented");
1454     }
1455    
1456    
1457     void RipleyDomain::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1458     {
1459     switch (functionSpaceType) {
1460     case Nodes:
1461     m_nodes->setTags(newTag, mask);
1462     break;
1463     case ReducedNodes:
1464     throw RipleyException("ReducedNodes does not support tags");
1465     case DegreesOfFreedom:
1466     throw RipleyException("DegreesOfFreedom does not support tags");
1467     case ReducedDegreesOfFreedom:
1468     throw RipleyException("ReducedDegreesOfFreedom does not support tags");
1469     case Elements:
1470     m_elements->setTags(newTag, mask);
1471     break;
1472     case ReducedElements:
1473     m_elements->setTags(newTag, mask);
1474     break;
1475     case FaceElements:
1476     m_faceElements->setTags(newTag, mask);
1477     break;
1478     case ReducedFaceElements:
1479     m_faceElements->setTags(newTag, mask);
1480     break;
1481     case Points:
1482     m_points->setTags(newTag, mask);
1483     break;
1484     default:
1485     stringstream msg;
1486     msg << "Ripley does not know anything about function space type " << functionSpaceType;
1487     throw RipleyException(msg.str());
1488     }
1489     }
1490    
1491     void RipleyDomain::setTagMap(const string& name, int tag)
1492     {
1493     m_tagMap[name] = tag;
1494     }
1495    
1496     int RipleyDomain::getTag(const string& name) const
1497     {
1498     return m_tagMap.find(name)->second;
1499     }
1500    
1501     bool RipleyDomain::isValidTagName(const string& name) const
1502     {
1503     return (m_tagMap.find(name)!=m_tagMap.end());
1504     }
1505    
1506     string RipleyDomain::showTagNames() const
1507     {
1508     stringstream ret;
1509     TagMap::const_iterator it;
1510     for (it=m_tagMap.begin(); it!=m_tagMap.end(); it++) {
1511     if (it!=m_tagMap.begin()) ret << ", ";
1512     ret << it->first;
1513     }
1514     return ret.str();
1515     }
1516    
1517     int RipleyDomain::getNumberOfTagsInUse(int functionSpaceCode) const
1518     {
1519     dim_t numTags=0;
1520     switch (functionSpaceCode) {
1521     case Nodes:
1522     numTags=m_nodes->getNumberOfTagsInUse();
1523     break;
1524     case ReducedNodes:
1525     throw RipleyException("ReducedNodes does not support tags");
1526     case DegreesOfFreedom:
1527     throw RipleyException("DegreesOfFreedom does not support tags");
1528     case ReducedDegreesOfFreedom:
1529     throw RipleyException("ReducedDegreesOfFreedom does not support tags");
1530     case Elements:
1531     case ReducedElements:
1532     numTags=m_elements->getNumberOfTagsInUse();
1533     break;
1534     case FaceElements:
1535     case ReducedFaceElements:
1536     numTags=m_faceElements->getNumberOfTagsInUse();
1537     break;
1538     case Points:
1539     numTags=m_points->getNumberOfTagsInUse();
1540     break;
1541     default:
1542     stringstream msg;
1543     msg << "Ripley does not know anything about function space type " << functionSpaceCode;
1544     throw RipleyException(msg.str());
1545     }
1546     return numTags;
1547     }
1548    
1549     const int* RipleyDomain::borrowListOfTagsInUse(int functionSpaceCode) const
1550     {
1551     switch (functionSpaceCode) {
1552     case Nodes:
1553     return &m_nodes->getTagsInUse()[0];
1554     case ReducedNodes:
1555     throw RipleyException("ReducedNodes does not support tags");
1556     case DegreesOfFreedom:
1557     throw RipleyException("DegreesOfFreedom does not support tags");
1558     case ReducedDegreesOfFreedom:
1559     throw RipleyException("ReducedDegreesOfFreedom does not support tags");
1560     case Elements:
1561     case ReducedElements:
1562     return &m_elements->getTagsInUse()[0];
1563     case FaceElements:
1564     case ReducedFaceElements:
1565     return &m_faceElements->getTagsInUse()[0];
1566     case Points:
1567     return &m_points->getTagsInUse()[0];
1568     default:
1569     stringstream msg;
1570     msg << "Ripley does not know anything about function space type " << functionSpaceCode;
1571     throw RipleyException(msg.str());
1572     }
1573     }
1574    
1575     bool RipleyDomain::canTag(int functionSpaceCode) const
1576     {
1577     switch(functionSpaceCode) {
1578     case Nodes:
1579     case Elements:
1580     case ReducedElements:
1581     case FaceElements:
1582     case ReducedFaceElements:
1583     case Points:
1584     return true;
1585     }
1586     return false;
1587     }
1588    
1589     escript::AbstractDomain::StatusType RipleyDomain::getStatus() const
1590     {
1591     return m_nodes->getStatus();
1592     }
1593    
1594     void RipleyDomain::prepare(bool optimize)
1595     {
1596     resolveNodeIds();
1597    
1598     // first step is to distribute the elements according to a global
1599     // distribution of DOF
1600     // - create dense labeling for the DOFs
1601     dim_t newGlobalNumDOFs=m_nodes->createDenseDOFLabeling();
1602    
1603     // - create a distribution of the global DOFs and determine
1604     // the MPI_rank controlling the DOFs on this processor
1605     IndexVector distribution(m_mpiInfo->size+1);
1606     Esys_MPIInfo_setDistribution(m_mpiInfo, 0, newGlobalNumDOFs-1, &distribution[0]);
1607     checkPasoError();
1608    
1609     // Now the mesh is redistributed according to the mpiRankOfDOF vector.
1610     // This will redistribute the Nodes and Elements including overlap and will
1611     // create an element coloring but will not create any mappings (see later
1612     // in this function).
1613     distributeByRankOfDOF(distribution);
1614    
1615     // At this stage we are able to start an optimization of the DOF
1616     // distribution using ParMetis. On return distribution is altered and new
1617     // DOF IDs have been assigned.
1618     if (optimize) {
1619     if (m_mpiInfo->size > 1) {
1620     optimizeDOFDistribution(distribution);
1621     distributeByRankOfDOF(distribution);
1622     }
1623    
1624     // optimize the local labeling of the degrees of freedom
1625     optimizeDOFLabeling(distribution);
1626     }
1627    
1628     // Rearrange elements in order to try and bring them closer to memory
1629     // locations of the nodes (distributed shared memory).
1630     optimizeElementOrdering();
1631    
1632     /* useful DEBUG:
1633     {
1634     printf("prepare: global DOF : %d\n",newGlobalNumDOFs);
1635     IndexPair range = m_nodes->getGlobalIdRange();
1636     printf("prepare: global node id range = %d :%d\n",range.first,range.second);
1637     range = m_nodes->getIdRange();
1638     printf("prepare: local node id range = %d :%d\n",range.first,range.second);
1639     }
1640     */
1641    
1642     // create the global indices
1643     IndexVector maskReducedNodes(m_nodes->getNumNodes(), -1);
1644     markNodes(maskReducedNodes, 0);
1645     IndexVector indexReducedNodes = packMask(maskReducedNodes);
1646    
1647     IndexVector nodeDistribution(m_mpiInfo->size+1);
1648     m_nodes->createDenseNodeLabeling(nodeDistribution, distribution);
1649     m_nodes->createDenseReducedDOFLabeling(maskReducedNodes);
1650     m_nodes->createDenseReducedNodeLabeling(maskReducedNodes);
1651    
1652     // create the missing mappings
1653     m_nodes->createNodeFileMappings(indexReducedNodes, distribution,
1654     nodeDistribution);
1655    
1656     updateTagsInUse();
1657     }
1658    
1659     void RipleyDomain::optimizeElementOrdering()
1660     {
1661     m_elements->optimizeOrdering();
1662     m_faceElements->optimizeOrdering();
1663     m_points->optimizeOrdering();
1664     }
1665    
1666     void RipleyDomain::updateTagsInUse()
1667     {
1668     m_nodes->updateTagsInUse();
1669     m_elements->updateTagsInUse();
1670     m_faceElements->updateTagsInUse();
1671     m_points->updateTagsInUse();
1672     }
1673    
1674     void RipleyDomain::createColoring(const IndexVector &node_localDOF_map)
1675     {
1676     m_elements->createColoring(node_localDOF_map);
1677     m_faceElements->createColoring(node_localDOF_map);
1678     m_points->createColoring(node_localDOF_map);
1679     }
1680    
1681     void RipleyDomain::distributeByRankOfDOF(const IndexVector &dofDistribution)
1682     {
1683     RankVector mpiRankOfDOF = m_nodes->getOwnerOfDOFs(dofDistribution);
1684    
1685     // First the elements are redistributed according to mpiRankOfDOF.
1686     // At the input the Node tables refer to the local labeling of
1687     // the nodes while at the output they refer to the global labeling
1688     // which is rectified in the next step
1689     m_elements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());
1690     m_faceElements->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());
1691     m_points->distributeByRankOfDOF(mpiRankOfDOF, m_nodes->getIdVector());
1692    
1693     // resolve the node IDs
1694     resolveNodeIds();
1695    
1696     // create a local labeling of the DOFs
1697     IndexPair dofRange = m_nodes->getDOFRange();
1698     dim_t len = dofRange.second - dofRange.first + 1;
1699    
1700     // local mask for used nodes
1701     IndexVector tmp_node_localDOF_mask(len, -1);
1702     IndexVector tmp_node_localDOF_map(m_nodes->getNumNodes(), -1);
1703    
1704     #pragma omp parallel for schedule(static)
1705     for (dim_t n = 0; n < m_nodes->getNumNodes(); n++) {
1706     #ifdef BOUNDS_CHECK
1707     if ((m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first) >= len
1708     || (m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first) < 0) {
1709     printf("BOUNDS_CHECK %s %d\n", __FILE__, __LINE__);
1710     exit(1);
1711     }
1712     #endif
1713     tmp_node_localDOF_mask[m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first] = n;
1714     }
1715    
1716     dim_t numDOFs = 0;
1717     for (dim_t n = 0; n < len; n++) {
1718     if (tmp_node_localDOF_mask[n] >= 0) {
1719     tmp_node_localDOF_mask[n] = numDOFs;
1720     numDOFs++;
1721     }
1722     }
1723     #pragma omp parallel for
1724     for (dim_t n = 0; n < m_nodes->getNumNodes(); n++) {
1725     tmp_node_localDOF_map[n] = tmp_node_localDOF_mask[m_nodes->getGlobalDegreesOfFreedom()[n] - dofRange.first];
1726     }
1727     // create element coloring
1728     createColoring(tmp_node_localDOF_map);
1729     }
1730    
1731     /**************************************************************
1732     Check whether there is any node which has no vertex. In case
1733     such a node exists, we don't use ParMetis since it requires
1734     that every node has at least 1 vertex (at line 129 of file
1735     "xyzpart.c" in parmetis 3.1.1, variable "nvtxs" would be 0 if
1736     any node has no vertex).
1737     **************************************************************/
1738     #ifdef USE_PARMETIS
1739     static int Check_Inputs_For_Parmetis(dim_t mpiSize, dim_t rank,
1740     const vector<dim_t> &distribution, MPI_Comm *comm)
1741     {
1742     dim_t i, len;
1743     int ret_val = 1;
1744    
1745     if (rank == 0) {
1746     for (i = 0; i < mpiSize; i++) {
1747     len = distribution[i + 1] - distribution[i];
1748     if (len == 0) {
1749     ret_val = 0;
1750     break;
1751     }
1752     }
1753     }
1754     MPI_Bcast(&ret_val, 1, MPI_INTEGER, 0, *comm);
1755     if (ret_val == 0)
1756     printf("INFO: Not using ParMetis since some nodes have no vertex!\n");
1757     return ret_val;
1758     }
1759     #endif
1760    
1761     void RipleyDomain::optimizeDOFDistribution(RankVector &distribution)
1762     {
1763     const Esys_MPI_rank myRank = m_mpiInfo->rank;
1764     dim_t mpiSize = m_mpiInfo->size;
1765     dim_t dim = getDim();
1766    
1767     // first step is to distribute the elements according to a global X of DOF
1768     const index_t myFirstVertex = distribution[myRank];
1769     const index_t myLastVertex = distribution[myRank + 1];
1770     const dim_t myNumVertices = myLastVertex - myFirstVertex;
1771     const dim_t globalNumVertices = distribution[mpiSize];
1772     dim_t len = 0;
1773     for (dim_t p = 0; p < mpiSize; ++p)
1774     len = MAX(len, distribution[p+1] - distribution[p]);
1775     index_t *partition = TMPMEMALLOC(len, index_t);
1776     float *xyz = TMPMEMALLOC(myNumVertices * dim, float);
1777     if (!(Esys_checkPtr(partition) || Esys_checkPtr(xyz))) {
1778     // set the coordinates
1779     // it is assumed that at least one node on this processor provides a
1780     // coordinate
1781     #pragma omp parallel for
1782     for (dim_t i = 0; i < m_nodes->getNumNodes(); ++i) {
1783     dim_t k = m_nodes->getGlobalDegreesOfFreedom()[i] - myFirstVertex;
1784     if ((k >= 0) && (k < myNumVertices)) {
1785     for (dim_t j = 0; j < dim; ++j)
1786     xyz[k * dim + j] = (float)(m_nodes->getCoordinates()[INDEX2(j, i, dim)]);
1787     }
1788     }
1789    
1790     IndexMatrix indexList(myNumVertices);
1791     // ksteube CSR of DOF IDs
1792     // create the adjacency structure xadj and adjncy
1793     // insert contributions from element matrices into columns of indexList
1794     m_elements->insertIntoIndexMatrixNoMainDiagonal(indexList,
1795     m_nodes->getGlobalDegreesOfFreedom(),
1796     m_nodes->getGlobalDegreesOfFreedom(),
1797     myFirstVertex, myLastVertex);
1798     m_faceElements->insertIntoIndexMatrixNoMainDiagonal(indexList,
1799     m_nodes->getGlobalDegreesOfFreedom(),
1800     m_nodes->getGlobalDegreesOfFreedom(),
1801     myFirstVertex, myLastVertex);
1802     m_points->insertIntoIndexMatrixNoMainDiagonal(indexList,
1803     m_nodes->getGlobalDegreesOfFreedom(),
1804     m_nodes->getGlobalDegreesOfFreedom(),
1805     myFirstVertex, myLastVertex);
1806    
1807     // create the local matrix pattern
1808     Paso_Pattern *pattern = createPatternFromIndexMatrix(
1809     indexList, 0, myNumVertices, 0, globalNumVertices, 0);
1810    
1811     #ifdef USE_PARMETIS
1812     if (mpiSize > 1 && Check_Inputs_For_Parmetis(mpiSize, myRank, distribution, &(m_mpiInfo->comm)) > 0) {
1813     int wgtflag = 0;
1814     int numflag = 0; /* pattern->ptr is C style: starting from 0 instead of 1 */
1815     int ncon = 1;
1816     int edgecut;
1817     int options[2];
1818     float *tpwgts = TMPMEMALLOC(ncon * mpiSize, float);
1819     float *ubvec = TMPMEMALLOC(ncon, float);
1820     for (dim_t i = 0; i < ncon * mpiSize; i++)
1821     tpwgts[i] = 1.0 / (float)mpiSize;
1822     for (dim_t i = 0; i < ncon; i++)
1823     ubvec[i] = 1.05;
1824     options[0] = 3;
1825     options[1] = 15;
1826     ParMETIS_V3_PartGeomKway(&distribution[0], pattern->ptr,
1827     pattern->index, NULL, NULL, &wgtflag, &numflag,
1828     &dim, xyz, &ncon, &mpiSize, tpwgts, ubvec, options,
1829     &edgecut, partition, // new CPU ownership of elements
1830     &m_mpiInfo->comm);
1831     TMPMEMFREE(ubvec);
1832     TMPMEMFREE(tpwgts);
1833     } else {
1834     for (dim_t i = 0; i < myNumVertices; ++i)
1835     partition[i] = 0; // CPU 0 owns it
1836     }
1837     #else
1838     for (dim_t i = 0; i < myNumVertices; ++i)
1839     partition[i] = myRank; // CPU myRank owns it
1840     #endif
1841    
1842     Paso_Pattern_free(pattern);
1843    
1844     // create a new distribution and labeling of the DOF
1845     RankVector newDistribution(mpiSize+1, 0);
1846     #pragma omp parallel
1847     {
1848     RankVector loc_partition_count(mpiSize, 0);
1849     #pragma omp for
1850     for (dim_t i = 0; i < myNumVertices; ++i)
1851     loc_partition_count[partition[i]]++;
1852     #pragma omp critical
1853     {
1854     for (dim_t i = 0; i < mpiSize; ++i)
1855     newDistribution[i] += loc_partition_count[i];
1856     }
1857     }
1858    
1859     // recvbuf will be the concatenation of each CPU's contribution
1860     // to newDistribution
1861     dim_t *recvbuf = TMPMEMALLOC(mpiSize * mpiSize, dim_t);
1862    
1863     #ifdef ESYS_MPI
1864     MPI_Allgather(&newDistribution[0], mpiSize, MPI_INT, recvbuf, mpiSize, MPI_INT, m_mpiInfo->comm);
1865     #else
1866     for (dim_t i = 0; i < mpiSize; ++i)
1867     recvbuf[i] = newDistribution[i];
1868     #endif
1869     newDistribution[0] = 0;
1870     IndexVector newGlobalDOFid(len);
1871     for (Esys_MPI_rank rank = 0; rank < mpiSize; rank++) {
1872     int c = 0;
1873     for (dim_t i = 0; i < myRank; ++i)
1874     c += recvbuf[rank + mpiSize * i];
1875     for (dim_t i = 0; i < myNumVertices; ++i) {
1876     if (rank == partition[i]) {
1877     newGlobalDOFid[i] = newDistribution[rank] + c;
1878     c++;
1879     }
1880     }
1881     for (dim_t i = myRank + 1; i < mpiSize; ++i)
1882     c += recvbuf[rank + mpiSize * i];
1883     newDistribution[rank + 1] = newDistribution[rank] + c;
1884     }
1885     TMPMEMFREE(recvbuf);
1886    
1887     // now the overlap needs to be created by sending the partition
1888     // around
1889     m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);
1890     for (dim_t i = 0; i < mpiSize+1; ++i)
1891     distribution[i] = newDistribution[i];
1892     }
1893     TMPMEMFREE(partition);
1894     TMPMEMFREE(xyz);
1895     }
1896    
1897     void RipleyDomain::optimizeDOFLabeling(const IndexVector &distribution)
1898     {
1899     Esys_MPI_rank myRank = m_mpiInfo->rank;
1900     index_t myFirstVertex = distribution[myRank];
1901     index_t myLastVertex = distribution[myRank + 1];
1902     dim_t myNumVertices = myLastVertex - myFirstVertex;
1903     dim_t len = 0;
1904     for (dim_t p = 0; p < m_mpiInfo->size; ++p)
1905     len = MAX(len, distribution[p + 1] - distribution[p]);
1906    
1907     IndexMatrix indexList(myNumVertices);
1908    
1909     // create the adjacency structure xadj and adjncy
1910     // insert contributions from element matrices into columns of indexList
1911     m_elements->insertIntoIndexMatrixNoMainDiagonal(indexList,
1912     m_nodes->getGlobalDegreesOfFreedom(),
1913     m_nodes->getGlobalDegreesOfFreedom(),
1914     myFirstVertex, myLastVertex);
1915     m_faceElements->insertIntoIndexMatrixNoMainDiagonal(indexList,
1916     m_nodes->getGlobalDegreesOfFreedom(),
1917     m_nodes->getGlobalDegreesOfFreedom(),
1918     myFirstVertex, myLastVertex);
1919     m_points->insertIntoIndexMatrixNoMainDiagonal(indexList,
1920     m_nodes->getGlobalDegreesOfFreedom(),
1921     m_nodes->getGlobalDegreesOfFreedom(),
1922     myFirstVertex, myLastVertex);
1923    
1924     // create the local matrix pattern
1925     Paso_Pattern *pattern = createPatternFromIndexMatrix(indexList, 0,
1926     myNumVertices, myFirstVertex, myLastVertex, -myFirstVertex);
1927    
1928     IndexVector newGlobalDOFid(len);
1929     Paso_Pattern_reduceBandwidth(pattern, &newGlobalDOFid[0]);
1930     Paso_Pattern_free(pattern);
1931    
1932     Esys_MPIInfo_noError(m_mpiInfo);
1933     checkPasoError();
1934    
1935     /* shift new labeling to create a global id */
1936     #pragma omp parallel for
1937     for (dim_t i = 0; i < myNumVertices; ++i)
1938     newGlobalDOFid[i] += myFirstVertex;
1939    
1940     /* distribute new labeling to other processors */
1941     m_nodes->resetGlobalDegreesOfFreedom(newGlobalDOFid, distribution);
1942     #if 0
1943     for (i = 0; i < m_nodes->getNumNodes(); ++i)
1944     printf("%d ", m_nodes->getGlobalDegreesOfFreedom()[i]);
1945     printf("\n");
1946     #endif
1947     }
1948    
1949     void RipleyDomain::resolveNodeIds()
1950     {
1951     // find the minimum and maximum node ID used by elements
1952     IndexPair minMaxId = m_elements->getNodeRange();
1953     index_t minId = minMaxId.first;
1954     index_t maxId = minMaxId.second;
1955     minMaxId = m_faceElements->getNodeRange();
1956     minId = MIN(minId, minMaxId.first);
1957     maxId = MAX(maxId, minMaxId.second);
1958     minMaxId = m_points->getNodeRange();
1959     minId = MIN(minId, minMaxId.first);
1960     maxId = MAX(maxId, minMaxId.second);
1961    
1962     #ifdef Ripley_TRACE
1963     index_t globalMinId, globalMaxId;
1964     #ifdef ESYS_MPI
1965     index_t idRange[2], globalIdRange[2];
1966     idRange[0] = -minId;
1967     idRange[1] = maxId;
1968     MPI_Allreduce(idRange, globalIdRange, 2, MPI_INT, MPI_MAX, m_mpiInfo->comm);
1969     globalMinId = -globalIdRange[0];
1970     globalMaxId = globalIdRange[1];
1971     #else
1972     globalMinId = minId;
1973     globalMaxId = maxId;
1974     #endif // ESYS_MPI
1975     cout << "Node id range used by elements is " << globalMinId << ":"
1976     << globalMaxId << endl;
1977     #endif // Ripley_TRACE
1978    
1979     // this is only true if we have no elements at all
1980     if (minId > maxId) {
1981     maxId = -1;
1982     minId = 0;
1983     }
1984    
1985     // allocate mappings for new local node labeling to global node labeling
1986     // (newLocalToGlobalNodeLabels) and global node labeling to the new local
1987     // node labeling. globalToNewLocalNodeLabels[i-minId] is the new local id
1988     // of global node i
1989     dim_t len = (maxId >= minId) ? maxId - minId + 1 : 0;
1990     IndexVector globalToNewLocalNodeLabels(len, -1);
1991    
1992     // mark the nodes referred by elements in globalToNewLocalNodeLabels
1993     // which is currently used as a mask
1994     markNodes(globalToNewLocalNodeLabels, minId);
1995    
1996     // create a local labeling newLocalToGlobalNodeLabels of the local
1997     // nodes by packing the mask globalToNewLocalNodeLabels
1998     IndexVector newLocalToGlobalNodeLabels = packMask(globalToNewLocalNodeLabels);
1999     const dim_t newNumNodes = newLocalToGlobalNodeLabels.size();
2000    
2001     // invert the new labeling and shift the index newLocalToGlobalNodeLabels
2002     // to global node IDs
2003     #pragma omp parallel for schedule(static)
2004     for (dim_t n = 0; n < newNumNodes; n++) {
2005     #ifdef BOUNDS_CHECK
2006     if (n >= len || n < 0) {
2007     printf("BOUNDS_CHECK %s %d n=%d\n", __FILE__, __LINE__, n);
2008     exit(1);
2009     }
2010     if (newLocalToGlobalNodeLabels[n] >= len || newLocalToGlobalNodeLabels[n] < 0) {
2011     printf("BOUNDS_CHECK %s %d n=%d\n", __FILE__, __LINE__, n);
2012     exit(1);
2013     }
2014     #endif
2015     globalToNewLocalNodeLabels[newLocalToGlobalNodeLabels[n]] = n;
2016     newLocalToGlobalNodeLabels[n] += minId;
2017     }
2018    
2019     // create a new node file
2020     m_nodes = m_nodes->gatherGlobal(newLocalToGlobalNodeLabels);
2021    
2022     // relabel nodes of the elements
2023     relabelElementNodes(globalToNewLocalNodeLabels, minId);
2024     }
2025    
2026     void RipleyDomain::relabelElementNodes(const IndexVector &newNode, index_t offset)
2027     {
2028     m_elements->relabelNodes(newNode, offset);
2029     m_faceElements->relabelNodes(newNode, offset);
2030     m_points->relabelNodes(newNode, offset);
2031     }
2032    
2033     void RipleyDomain::markNodes(IndexVector &mask, index_t offset)
2034     {
2035     m_elements->markNodes(mask, offset);
2036     m_faceElements->markNodes(mask, offset);
2037     m_points->markNodes(mask, offset);
2038     }
2039    
2040     void RipleyDomain::createMappings(const IndexVector &dofDistribution, const IndexVector &nodeDistribution)
2041     {
2042     IndexVector maskReducedNodes(m_nodes->getNumNodes(), -1);
2043     markNodes(maskReducedNodes, 0);
2044    
2045     IndexVector indexReducedNodes = packMask(maskReducedNodes);
2046     m_nodes->createNodeFileMappings(indexReducedNodes, dofDistribution,
2047     nodeDistribution);
2048     }
2049    
2050     Paso_SystemMatrixPattern *RipleyDomain::makePattern(bool reduce_row_order, bool reduce_col_order) const
2051     {
2052     Paso_Connector *colConnector, *rowConnector;
2053     NodeMapping_ptr colMap, rowMap;
2054     Paso_Distribution *colDistribution = NULL, *rowDistribution = NULL;
2055    
2056     if (reduce_col_order) {
2057     colMap = m_nodes->getReducedDegreesOfFreedomMapping();
2058     colDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();
2059     colConnector = m_nodes->getReducedDegreesOfFreedomConnector();
2060     } else {
2061     colMap = m_nodes->getDegreesOfFreedomMapping();
2062     colDistribution = m_nodes->getDegreesOfFreedomDistribution();
2063     colConnector = m_nodes->getDegreesOfFreedomConnector();
2064     }
2065    
2066     if (reduce_row_order) {
2067     rowMap = m_nodes->getReducedDegreesOfFreedomMapping();
2068     rowDistribution = m_nodes->getReducedDegreesOfFreedomDistribution();
2069     rowConnector = m_nodes->getReducedDegreesOfFreedomConnector();
2070     } else {
2071     rowMap = m_nodes->getDegreesOfFreedomMapping();
2072     rowDistribution = m_nodes->getDegreesOfFreedomDistribution();
2073     rowConnector = m_nodes->getDegreesOfFreedomConnector();
2074     }
2075    
2076     // insert contributions from element matrices into columns in indexList
2077     IndexMatrix indexList(rowMap->numTargets);
2078     m_elements->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);
2079     m_faceElements->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);
2080     m_points->insertIntoIndexMatrix(indexList, rowMap->target, colMap->target);
2081    
2082     // create patterns
2083     Paso_Pattern *mainPattern = createPatternFromIndexMatrix(indexList, 0,
2084     Paso_Distribution_getMyNumComponents(rowDistribution),
2085     0, Paso_Distribution_getMyNumComponents(colDistribution), 0);
2086     Paso_Pattern *colCouplePattern = createPatternFromIndexMatrix(indexList, 0,
2087     Paso_Distribution_getMyNumComponents(rowDistribution),
2088     Paso_Distribution_getMyNumComponents(colDistribution),
2089     colMap->numTargets,
2090     -Paso_Distribution_getMyNumComponents(colDistribution));
2091     Paso_Pattern *rowCouplePattern = createPatternFromIndexMatrix(indexList,
2092     Paso_Distribution_getMyNumComponents(rowDistribution),
2093     rowMap->numTargets, 0,
2094     Paso_Distribution_getMyNumComponents(colDistribution), 0);
2095    
2096     // if everything is in order we can create the return value
2097     Paso_SystemMatrixPattern *out = Paso_SystemMatrixPattern_alloc(
2098     MATRIX_FORMAT_DEFAULT, rowDistribution, colDistribution,
2099     mainPattern, colCouplePattern, rowCouplePattern,
2100     colConnector, rowConnector);
2101    
2102     // clean up
2103     Paso_Pattern_free(mainPattern);
2104     Paso_Pattern_free(colCouplePattern);
2105     Paso_Pattern_free(rowCouplePattern);
2106     Esys_MPIInfo_noError(m_mpiInfo);
2107     return out;
2108     }
2109    
2110     Paso_SystemMatrixPattern *RipleyDomain::getPattern(bool reduce_row_order, bool reduce_col_order) const
2111     {
2112     return Paso_SystemMatrixPattern_getReference(makePattern(reduce_row_order, reduce_col_order));
2113     /*
2114     Paso_SystemMatrixPattern *out = NULL;
2115    
2116     // make sure that the requested pattern is available
2117     if (reduce_row_order) {
2118     if (reduce_col_order) {
2119     if (m_reducedReducedPattern == NULL)
2120     m_reducedReducedPattern = makePattern(reduce_row_order, reduce_col_order);
2121     } else {
2122     if (m_reducedFullPattern == NULL)
2123     m_reducedFullPattern = makePattern(reduce_row_order, reduce_col_order);
2124     }
2125     } else {
2126     if (reduce_col_order) {
2127     if (m_fullReducedPattern == NULL)
2128     m_fullReducedPattern = makePattern(reduce_row_order, reduce_col_order);
2129     } else {
2130     if (m_fullFullPattern == NULL)
2131     m_fullFullPattern = makePattern(reduce_row_order, reduce_col_order);
2132     }
2133     }
2134     if (Ripley_noError()) {
2135     if (reduce_row_order) {
2136     if (reduce_col_order) {
2137     out = Paso_SystemMatrixPattern_getReference(m_reducedReducedPattern);
2138     } else {
2139     out = Paso_SystemMatrixPattern_getReference(m_reducedFullPattern);
2140     }
2141     } else {
2142     if (reduce_col_order) {
2143     out = Paso_SystemMatrixPattern_getReference(m_fullReducedPattern);
2144     } else {
2145     out = Paso_SystemMatrixPattern_getReference(m_fullFullPattern);
2146     }
2147     }
2148     }
2149     return out;
2150     */
2151     }
2152    
2153     } // end of namespace ripley
2154    

  ViewVC Help
Powered by ViewVC 1.1.26