00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
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 : setupVerb_(setupVerb),
00116 map_(),
00117 mesh_(mesh),
00118 subdomains_(),
00119 basis_(),
00120 vecSpace_(),
00121 vecType_(vecType),
00122 ghostImporter_()
00123 ,transformationBuilder_(0)
00124 {
00125 init(tuple(funcDomains), BasisArray(tuple(basis)));
00126 }
00127
00128
00129 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00130 const CellFilter& funcDomains,
00131 const VectorType<double>& vecType,
00132 int setupVerb)
00133 : setupVerb_(setupVerb),
00134 map_(),
00135 mesh_(mesh),
00136 subdomains_(),
00137 basis_(),
00138 vecSpace_(),
00139 vecType_(vecType),
00140 ghostImporter_()
00141 ,transformationBuilder_(0)
00142 {
00143 init(Array<CellFilter>(basis.size(), funcDomains), basis);
00144 }
00145
00146
00147
00148 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00149 const RCP<DOFMapBase>& map,
00150 const RCP<Array<int> >& bcIndices,
00151 const VectorType<double>& vecType,
00152 int setupVerb)
00153 : setupVerb_(setupVerb),
00154 map_(map),
00155 mesh_(mesh),
00156 subdomains_(),
00157 basis_(),
00158 vecSpace_(),
00159 vecType_(vecType),
00160 ghostImporter_()
00161 ,transformationBuilder_(0)
00162 {
00163 init(map->funcDomains(), basis, bcIndices, true);
00164 }
00165
00166 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00167 const RCP<DOFMapBase>& map,
00168 const VectorType<double>& vecType,
00169 int setupVerb)
00170 : map_(map),
00171 mesh_(mesh),
00172 subdomains_(),
00173 basis_(),
00174 vecSpace_(),
00175 vecType_(vecType),
00176 ghostImporter_()
00177 ,transformationBuilder_(0)
00178 {
00179 init(map->funcDomains(), basis);
00180 }
00181
00182
00183 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisFamily& basis,
00184 const SpectralBasis& spBasis,
00185 const VectorType<double>& vecType,
00186 int setupVerb)
00187 : setupVerb_(setupVerb),
00188 map_(),
00189 mesh_(mesh),
00190 subdomains_(),
00191 basis_(),
00192 vecSpace_(),
00193 vecType_(vecType),
00194 ghostImporter_()
00195 ,transformationBuilder_(0)
00196 {
00197 init(maximalRegions(spBasis.nterms()),
00198 replicate(basis, spBasis.nterms()));
00199 }
00200
00201 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00202 const SpectralBasis& spBasis,
00203 const VectorType<double>& vecType,
00204 int setupVerb)
00205 : setupVerb_(setupVerb),
00206 map_(),
00207 mesh_(mesh),
00208 subdomains_(),
00209 basis_(),
00210 vecSpace_(),
00211 vecType_(vecType),
00212 ghostImporter_()
00213 ,transformationBuilder_(0)
00214 {
00215 init(maximalRegions(basis.size() * spBasis.nterms()),
00216 replicate(basis, spBasis.nterms()));
00217 }
00218
00219 DiscreteSpace::DiscreteSpace(const Mesh& mesh, const BasisArray& basis,
00220 const RCP<FunctionSupportResolver>& fsr,
00221 const VectorType<double>& vecType,
00222 int setupVerb)
00223 : setupVerb_(setupVerb),
00224 map_(),
00225 mesh_(mesh),
00226 subdomains_(),
00227 basis_(basis),
00228 vecSpace_(),
00229 vecType_(vecType),
00230 ghostImporter_()
00231 ,transformationBuilder_(new DiscreteSpaceTransfBuilder())
00232 {
00233 bool partitionBCs = false;
00234 DOFMapBuilder builder(mesh, fsr, partitionBCs, setupVerb);
00235
00236 map_ = builder.colMap()[0];
00237 Array<Set<CellFilter> > cf = builder.unkCellFilters()[0];
00238
00239 for (int i=0; i<cf.size(); i++)
00240 {
00241 Array<Array<CellFilter> > dimCF(mesh.spatialDim()+1);
00242 for (Set<CellFilter>::const_iterator
00243 iter=cf[i].begin(); iter != cf[i].end(); iter++)
00244 {
00245 CellFilter f = *iter;
00246 int dim = f.dimension(mesh);
00247 dimCF[dim].append(f);
00248 }
00249 for (int d=mesh.spatialDim(); d>=0; d--)
00250 {
00251 if (dimCF[d].size() == 0) continue;
00252 CellFilter f = dimCF[d][0];
00253 for (int j=1; j<dimCF[d].size(); j++)
00254 {
00255 f = f + dimCF[d][j];
00256 }
00257 subdomains_.append(f);
00258 break;
00259 }
00260 }
00261 RCP<Array<int> > dummyBCIndices;
00262
00263
00264 transformationBuilder_ = rcp(new DiscreteSpaceTransfBuilder( mesh , basis , map_ ));
00265
00266 initVectorSpace(dummyBCIndices, partitionBCs);
00267 initImporter();
00268 }
00269
00270 void DiscreteSpace::init(
00271 const Array<CellFilter>& regions,
00272 const BasisArray& basis)
00273 {
00274 RCP<Array<int> > dummyBCIndices;
00275 init(regions, basis, dummyBCIndices, false);
00276 }
00277
00278 void DiscreteSpace::init(
00279 const Array<CellFilter>& regions,
00280 const BasisArray& basis,
00281 const RCP<Array<int> >& isBCIndex,
00282 bool partitionBCs)
00283 {
00284 basis_ = basis;
00285 subdomains_ = regions;
00286 Array<RCP<BasisDOFTopologyBase> > basisTop(basis.size());
00287 for (int b=0; b<basis.size(); b++)
00288 {
00289 basisTop[b] = rcp_dynamic_cast<BasisDOFTopologyBase>(basis[b].ptr());
00290 }
00291
00292 if (map_.get()==0)
00293 {
00294 Array<Set<CellFilter> > cf(regions.size());
00295 for (int i=0; i<regions.size(); i++) cf[i] = makeSet(regions[i]);
00296 DOFMapBuilder b(setupVerb_);
00297 map_ = b.makeMap(mesh_, basisTop, cf);
00298 }
00299
00300
00301 transformationBuilder_ = rcp(new DiscreteSpaceTransfBuilder( mesh_ , basis , map_ ));
00302
00303 initVectorSpace(isBCIndex, partitionBCs);
00304
00305 if (!partitionBCs)
00306 {
00307 initImporter();
00308 }
00309 }
00310
00311 void DiscreteSpace::initVectorSpace(
00312 const RCP<Array<int> >& isBCIndex,
00313 bool partitionBCs)
00314 {
00315 TEUCHOS_TEST_FOR_EXCEPTION(map_.get()==0, std::logic_error,
00316 "uninitialized map");
00317
00318 int nDof = map_->numLocalDOFs();
00319 int lowDof = map_->lowestLocalDOF();
00320
00321 if (partitionBCs)
00322 {
00323 TEUCHOS_TEST_FOR_EXCEPT(isBCIndex.get() == 0);
00324
00325 int nBCDofs = 0;
00326 for (int i=0; i<nDof; i++)
00327 {
00328 if ((*isBCIndex)[i]) nBCDofs++;
00329 }
00330
00331 int nTotalBCDofs = nBCDofs;
00332 mesh().comm().allReduce(&nBCDofs, &nTotalBCDofs, 1, MPIDataType::intType(), MPIOp::sumOp());
00333 int nTotalInteriorDofs = map_->numDOFs() - nTotalBCDofs;
00334
00335 Array<int> interiorDofs(nDof - nBCDofs);
00336 Array<int> bcDofs(nBCDofs);
00337 int iBC = 0;
00338 int iIn = 0;
00339 for (int i=0; i<nDof; i++)
00340 {
00341 if ((*isBCIndex)[i]) bcDofs[iBC++] = lowDof+i;
00342 else interiorDofs[iIn++] = lowDof+i;
00343 }
00344 const int* bcDofPtr = vecPtr(bcDofs);
00345 const int* intDofPtr = vecPtr(interiorDofs);
00346 VectorSpace<double> bcSpace = vecType_.createSpace(nTotalBCDofs, nBCDofs,
00347 bcDofPtr, mesh().comm());
00348 VectorSpace<double> interiorSpace = vecType_.createSpace(nTotalInteriorDofs, nDof-nBCDofs,
00349 intDofPtr, mesh().comm());
00350
00351 vecSpace_ = blockSpace<double>(interiorSpace, bcSpace);
00352 }
00353 else
00354 {
00355 Array<int> dofs(nDof);
00356 for (int i=0; i<nDof; i++) dofs[i] = lowDof + i;
00357
00358 vecSpace_ = vecType_.createSpace(map_->numDOFs(),
00359 map_->numLocalDOFs(),
00360 &(dofs[0]), mesh().comm());
00361 }
00362 }
00363
00364 void DiscreteSpace::initImporter()
00365 {
00366 TEUCHOS_TEST_FOR_EXCEPTION(map_.get()==0, std::logic_error,
00367 "uninitialized map");
00368 TEUCHOS_TEST_FOR_EXCEPTION(vecSpace_.ptr().get()==0, std::logic_error,
00369 "uninitialized vector space");
00370 TEUCHOS_TEST_FOR_EXCEPTION(vecType_.ptr().get()==0, std::logic_error,
00371 "uninitialized vector type");
00372
00373 RCP<Array<int> > ghostIndices = map_->ghostIndices();
00374 int nGhost = ghostIndices->size();
00375 int* ghosts = 0;
00376 if (nGhost!=0) ghosts = &((*ghostIndices)[0]);
00377 ghostImporter_ = vecType_.createGhostImporter(vecSpace_, nGhost, ghosts);
00378 }
00379
00380
00381 Array<CellFilter> DiscreteSpace::maximalRegions(int n) const
00382 {
00383 CellFilter cf = new MaximalCellFilter();
00384 Array<CellFilter> rtn(n, cf);
00385 return rtn;
00386 }
00387
00388
00389
00390 void DiscreteSpace::importGhosts(const Vector<double>& x,
00391 RCP<GhostView<double> >& ghostView) const
00392 {
00393 ghostImporter_->importView(x, ghostView);
00394 }