Sundance 2.1

00001 /* @HEADER@ */
00002 // ************************************************************************
00003 // 
00004 //                              Sundance
00005 //                 Copyright (2005) Sandia Corporation
00006 // 
00007 // Copyright (year first published) Sandia Corporation.  Under the terms 
00008 // of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 
00009 // retains certain rights in this software.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //                                                                                 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA                                                                                
00025 // Questions? Contact Kevin Long (krlong@sandia.gov), 
00026 // Sandia National Laboratories, Livermore, California, USA
00027 // 
00028 // ************************************************************************
00029 /* @HEADER@ */
00030 
00031 #include "SundanceDiscreteSpace.hpp"
00032 #include "SundanceDOFMapBuilder.hpp"
00033 #include "SundanceFunctionSupportResolver.hpp"
00034 #include "SundanceOut.hpp"
00035 #include "SundanceMaximalCellFilter.hpp"
00036 #include "PlayaDefaultBlockVectorSpaceDecl.hpp"
00037 #include "PlayaMPIContainerComm.hpp"
00038 
00039 
00040 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00041 #include "PlayaVectorImpl.hpp"
00042 #include "PlayaDefaultBlockVectorImpl.hpp"
00043 #endif
00044 
00045 
00046 using namespace Sundance;
00047 using namespace Teuchos;
00048 using Playa::blockSpace;
00049 using Playa::MPIComm;
00050 using Playa::MPIContainerComm;
00051 using Playa::MPIDataType;
00052 using Playa::MPIOp;
00053 
00054 const int* vecPtr(const Array<int>& x)
00055 {
00056   static const int* dum = 0;
00057   if (x.size()==0) return dum;
00058   else return &(x[0]);
00059 }
00060 
00061 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisFamily& basis,
00062   const VectorType<double>& vecType,
00063   int setupVerb)
00064   : setupVerb_(setupVerb), 
00065     map_(),
00066     mesh_(mesh), 
00067     subdomains_(),
00068     basis_(),
00069     vecSpace_(), 
00070     vecType_(vecType),
00071     ghostImporter_()
00072   ,transformationBuilder_(0)
00073 {
00074   init(maximalRegions(1), BasisArray(tuple(basis)));
00075 }
00076 
00077 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00078   const VectorType<double>& vecType,
00079   int setupVerb)
00080   : setupVerb_(setupVerb),
00081     map_(), 
00082     mesh_(mesh), 
00083     subdomains_(),
00084     basis_(),
00085     vecSpace_(), 
00086     vecType_(vecType),
00087     ghostImporter_()
00088   ,transformationBuilder_(0)
00089 {
00090   init(maximalRegions(basis.size()), basis);
00091 }
00092 
00093 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00094   const Array<CellFilter>& funcDomains,
00095   const VectorType<double>& vecType,
00096   int setupVerb)
00097   : setupVerb_(setupVerb),
00098     map_(), 
00099     mesh_(mesh), 
00100     subdomains_(),
00101     basis_(),
00102     vecSpace_(), 
00103     vecType_(vecType),
00104     ghostImporter_()
00105   ,transformationBuilder_(0)
00106 {
00107   init(funcDomains, basis);
00108 }
00109 
00110 
00111 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisFamily& basis,
00112   const CellFilter& funcDomains,
00113   const VectorType<double>& vecType,
00114   int setupVerb)
00115   : map_(), 
00116     mesh_(mesh), 
00117     subdomains_(),
00118     basis_(),
00119     vecSpace_(), 
00120     vecType_(vecType),
00121     ghostImporter_()
00122   ,transformationBuilder_(0)
00123 {
00124   init(tuple(funcDomains), BasisArray(tuple(basis)));
00125 }
00126 
00127 
00128 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00129   const CellFilter& funcDomains,
00130   const VectorType<double>& vecType,
00131   int setupVerb)
00132   : setupVerb_(setupVerb),
00133     map_(), 
00134     mesh_(mesh), 
00135     subdomains_(),
00136     basis_(),
00137     vecSpace_(), 
00138     vecType_(vecType),
00139     ghostImporter_()
00140   ,transformationBuilder_(0)
00141 {
00142   init(Array<CellFilter>(basis.size(), funcDomains), basis);
00143 }
00144 
00145 
00146 
00147 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00148   const RCP<DOFMapBase>& map,
00149   const RCP<Array<int> >& bcIndices,
00150   const VectorType<double>& vecType,
00151   int setupVerb)
00152   : setupVerb_(setupVerb),
00153     map_(map), 
00154     mesh_(mesh),
00155     subdomains_(),
00156     basis_(),
00157     vecSpace_(), 
00158     vecType_(vecType),
00159     ghostImporter_()
00160   ,transformationBuilder_(0)
00161 {
00162   init(map->funcDomains(), basis, bcIndices, true);
00163 }
00164 
00165 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00166   const RCP<DOFMapBase>& map,
00167   const VectorType<double>& vecType,
00168   int setupVerb)
00169   : map_(map), 
00170     mesh_(mesh),
00171     subdomains_(),
00172     basis_(),
00173     vecSpace_(), 
00174     vecType_(vecType),
00175     ghostImporter_()
00176   ,transformationBuilder_(0)
00177 {
00178   init(map->funcDomains(), basis);
00179 }
00180 
00181 
00182 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisFamily& basis,
00183   const SpectralBasis& spBasis,
00184   const VectorType<double>& vecType,
00185   int setupVerb)
00186   : setupVerb_(setupVerb),
00187     map_(),
00188     mesh_(mesh), 
00189     subdomains_(),
00190     basis_(),
00191     vecSpace_(), 
00192     vecType_(vecType),
00193     ghostImporter_()
00194   ,transformationBuilder_(0)
00195 {
00196   init(maximalRegions(spBasis.nterms()), 
00197     replicate(basis, spBasis.nterms()));
00198 }
00199 
00200 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00201   const SpectralBasis& spBasis,
00202   const VectorType<double>& vecType,
00203   int setupVerb)
00204   : setupVerb_(setupVerb),
00205     map_(), 
00206     mesh_(mesh), 
00207     subdomains_(),
00208     basis_(),
00209     vecSpace_(), 
00210     vecType_(vecType),
00211     ghostImporter_()
00212   ,transformationBuilder_(0)
00213 {
00214   init(maximalRegions(basis.size() * spBasis.nterms()), 
00215     replicate(basis, spBasis.nterms()));
00216 }
00217 
00218 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00219   const RCP<FunctionSupportResolver>& fsr,
00220   const VectorType<double>& vecType,
00221   int setupVerb)
00222   : setupVerb_(setupVerb),
00223     map_(), 
00224     mesh_(mesh), 
00225     subdomains_(),
00226     basis_(basis),
00227     vecSpace_(), 
00228     vecType_(vecType),
00229     ghostImporter_()
00230   ,transformationBuilder_(new DiscreteSpaceTransfBuilder())
00231 {
00232   bool partitionBCs = false;
00233   DOFMapBuilder builder(mesh, fsr, partitionBCs, setupVerb);
00234 
00235   map_ = builder.colMap()[0];
00236   Array<Set<CellFilter> > cf = builder.unkCellFilters()[0];
00237 
00238   for (int i=0; i<cf.size(); i++)
00239   {
00240     Array<Array<CellFilter> > dimCF(mesh.spatialDim()+1);
00241     for (Set<CellFilter>::const_iterator 
00242            iter=cf[i].begin(); iter != cf[i].end(); iter++)
00243     {
00244       CellFilter f = *iter;
00245       int dim = f.dimension(mesh);
00246       dimCF[dim].append(f);
00247     }
00248     for (int d=mesh.spatialDim(); d>=0; d--)
00249     {
00250       if (dimCF[d].size() == 0) continue;
00251       CellFilter f = dimCF[d][0];
00252       for (int j=1; j<dimCF[d].size(); j++)
00253       {
00254         f = f + dimCF[d][j];
00255       }
00256       subdomains_.append(f);
00257       break;
00258     }
00259   }
00260   RCP<Array<int> > dummyBCIndices;
00261   
00262   // set up the transformation
00263   transformationBuilder_ = rcp(new DiscreteSpaceTransfBuilder( mesh , basis , map_ ));
00264 
00265   initVectorSpace(dummyBCIndices, partitionBCs);
00266   initImporter();
00267 }
00268 
00269 void DiscreteSpace::init(
00270   const Array<CellFilter>& regions,
00271   const BasisArray& basis)
00272 {
00273   RCP<Array<int> > dummyBCIndices;
00274   init(regions, basis, dummyBCIndices, false);
00275 }
00276 
00277 void DiscreteSpace::init(
00278   const Array<CellFilter>& regions,
00279   const BasisArray& basis,
00280   const RCP<Array<int> >& isBCIndex, 
00281   bool partitionBCs)
00282 {
00283   basis_ = basis;
00284   subdomains_ = regions;
00285   Array<RCP<BasisDOFTopologyBase> > basisTop(basis.size());
00286   for (int b=0; b<basis.size(); b++)
00287   {
00288     basisTop[b] = rcp_dynamic_cast<BasisDOFTopologyBase>(basis[b].ptr());
00289   }
00290 
00291   if (map_.get()==0) 
00292   {
00293     Array<Set<CellFilter> > cf(regions.size());
00294     for (int i=0; i<regions.size(); i++) cf[i] = makeSet(regions[i]);
00295     DOFMapBuilder b(setupVerb_);
00296     map_ = b.makeMap(mesh_, basisTop, cf);
00297   }
00298 
00299   // set up the transformation
00300   transformationBuilder_ = rcp(new DiscreteSpaceTransfBuilder( mesh_ , basis , map_ ));
00301 
00302   initVectorSpace(isBCIndex, partitionBCs);
00303 
00304   if (!partitionBCs)
00305   {
00306     initImporter();
00307   }
00308 }
00309 
00310 void DiscreteSpace::initVectorSpace(
00311   const RCP<Array<int> >& isBCIndex, 
00312   bool partitionBCs)
00313 {
00314   TEUCHOS_TEST_FOR_EXCEPTION(map_.get()==0, std::logic_error,
00315     "uninitialized map");
00316 
00317   int nDof = map_->numLocalDOFs();
00318   int lowDof = map_->lowestLocalDOF();
00319 
00320   if (partitionBCs)
00321   {
00322     TEUCHOS_TEST_FOR_EXCEPT(isBCIndex.get() == 0);
00323 
00324     int nBCDofs = 0;
00325     for (int i=0; i<nDof; i++)
00326     {
00327       if ((*isBCIndex)[i]) nBCDofs++;
00328     }
00329     
00330     int nTotalBCDofs = nBCDofs;
00331     mesh().comm().allReduce(&nBCDofs, &nTotalBCDofs, 1, MPIDataType::intType(), MPIOp::sumOp());
00332     int nTotalInteriorDofs = map_->numDOFs() - nTotalBCDofs;
00333 
00334     Array<int> interiorDofs(nDof - nBCDofs);
00335     Array<int> bcDofs(nBCDofs);
00336     int iBC = 0;
00337     int iIn = 0;
00338     for (int i=0; i<nDof; i++)
00339     {
00340       if ((*isBCIndex)[i]) bcDofs[iBC++] = lowDof+i;
00341       else interiorDofs[iIn++] = lowDof+i;
00342     }
00343     const int* bcDofPtr = vecPtr(bcDofs);
00344     const int* intDofPtr = vecPtr(interiorDofs);
00345     VectorSpace<double> bcSpace = vecType_.createSpace(nTotalBCDofs, nBCDofs,
00346       bcDofPtr, mesh().comm());
00347     VectorSpace<double> interiorSpace = vecType_.createSpace(nTotalInteriorDofs, nDof-nBCDofs,
00348       intDofPtr, mesh().comm());
00349 
00350     vecSpace_ = blockSpace<double>(interiorSpace, bcSpace);
00351   }
00352   else
00353   {
00354     Array<int> dofs(nDof);
00355     for (int i=0; i<nDof; i++) dofs[i] = lowDof + i;
00356     
00357     vecSpace_ = vecType_.createSpace(map_->numDOFs(),
00358       map_->numLocalDOFs(),
00359       &(dofs[0]), mesh().comm());
00360   }
00361 }
00362 
00363 void DiscreteSpace::initImporter()
00364 {
00365   TEUCHOS_TEST_FOR_EXCEPTION(map_.get()==0, std::logic_error,
00366     "uninitialized map");
00367   TEUCHOS_TEST_FOR_EXCEPTION(vecSpace_.ptr().get()==0, std::logic_error,
00368     "uninitialized vector space");
00369   TEUCHOS_TEST_FOR_EXCEPTION(vecType_.ptr().get()==0, std::logic_error,
00370     "uninitialized vector type");
00371   
00372   RCP<Array<int> > ghostIndices = map_->ghostIndices();
00373   int nGhost = ghostIndices->size();
00374   int* ghosts = 0;
00375   if (nGhost!=0) ghosts = &((*ghostIndices)[0]);
00376   ghostImporter_ = vecType_.createGhostImporter(vecSpace_, nGhost, ghosts);
00377 }
00378 
00379 
00380 Array<CellFilter> DiscreteSpace::maximalRegions(int n) const
00381 {
00382   CellFilter cf = new MaximalCellFilter();
00383   Array<CellFilter> rtn(n, cf);
00384   return rtn;
00385 }
00386 
00387 
00388 
00389 void DiscreteSpace::importGhosts(const Vector<double>& x,
00390   RCP<GhostView<double> >& ghostView) const
00391 {
00392   ghostImporter_->importView(x, ghostView);
00393 }
SundanceDiscreteSpace.cpp
