MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_BlockedCoordinatesTransferFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46#ifndef MUELU_BLOCKEDCOORDINATESTRANSFER_FACTORY_DEF_HPP
47#define MUELU_BLOCKEDCOORDINATESTRANSFER_FACTORY_DEF_HPP
48
49#include "Xpetra_ImportFactory.hpp"
50#include "Xpetra_MultiVectorFactory.hpp"
51#include "Xpetra_MapFactory.hpp"
52#include "Xpetra_IO.hpp"
53
54#include "MueLu_CoarseMapFactory.hpp"
55#include "MueLu_Aggregates.hpp"
57
58#include "MueLu_Level.hpp"
59#include "MueLu_Monitor.hpp"
60
61namespace MueLu {
62
63 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65 RCP<ParameterList> validParamList = rcp(new ParameterList());
66
67 validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Factory for coordinates generation");
68 validParamList->set<RCP<const FactoryBase> >("CoarseMap", Teuchos::null, "Generating factory of the coarse map");
69 return validParamList;
70 }
71
72 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
74 Input(coarseLevel, "CoarseMap");
75
76 // Make sure the Level knows I need these sub-Factories
77 const size_t numSubFactories = NumFactories();
78 for(size_t i=0; i<numSubFactories; i++) {
79 const RCP<const FactoryBase>& myFactory = subFactories_[i];
80 coarseLevel.DeclareInput("Coordinates", myFactory.getRawPtr(), this);
81 }
82
83 // call DeclareInput of all user-given transfer factories
84 for (std::vector<RCP<const FactoryBase> >::const_iterator it = subFactories_.begin(); it != subFactories_.end(); ++it)
85 (*it)->CallDeclareInput(coarseLevel);
86 }
87
88 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
90 FactoryMonitor m(*this, "Build", coarseLevel);
91
92 typedef Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> dMV;
93 typedef Xpetra::BlockedMultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> dBV;
94
95 GetOStream(Runtime0) << "Transferring (blocked) coordinates" << std::endl;
96
97 const size_t numSubFactories = NumFactories();
98 std::vector<RCP<const Map> > subBlockMaps(numSubFactories);
99 std::vector<RCP<dMV> > subBlockCoords(numSubFactories);
100
101 if (coarseLevel.IsAvailable("Coordinates", this)) {
102 GetOStream(Runtime0) << "Reusing coordinates" << std::endl;
103 return;
104 }
105
106 // Get components
107 for(size_t i=0; i<numSubFactories; i++) {
108 GetOStream(Runtime1) << "Generating Coordinates for block " << i <<"/"<<numSubFactories <<std::endl;
109 const RCP<const FactoryBase>& myFactory = subFactories_[i];
110 myFactory->CallBuild(coarseLevel);
111 subBlockCoords[i] = coarseLevel.Get<RCP<dMV> >("Coordinates", myFactory.get());
112 subBlockMaps[i] = subBlockCoords[i]->getMap();
113 }
114
115 // Blocked Map
116 RCP<const BlockedMap> coarseCoordMapBlocked;
117
118 {
119 // coarseMap is being used to set up the domain map of tentative P, and therefore, the row map of Ac
120 // Therefore, if we amalgamate coarseMap, logical nodes in the coordinates vector would correspond to
121 // logical blocks in the matrix
122 RCP<const BlockedMap> coarseMap = Get< RCP<const BlockedMap> >(coarseLevel, "CoarseMap");
123 bool thyraMode = coarseMap->getThyraMode();
124
125 ArrayView<const GO> elementAList = coarseMap->getFullMap()->getLocalElementList();
126
127 LO blkSize = 1;
128 if (rcp_dynamic_cast<const StridedMap>(coarseMap->getMap(0, thyraMode)) != Teuchos::null)
129 blkSize = rcp_dynamic_cast<const StridedMap>(coarseMap->getMap(0, thyraMode))->getFixedBlockSize();
130
131 for(size_t i=1; i<numSubFactories; i++) {
132 LO otherBlkSize = 1;
133 if (rcp_dynamic_cast<const StridedMap>(coarseMap->getMap(i, thyraMode)) != Teuchos::null)
134 otherBlkSize = rcp_dynamic_cast<const StridedMap>(coarseMap->getMap(i, thyraMode))->getFixedBlockSize();
135 TEUCHOS_TEST_FOR_EXCEPTION(otherBlkSize != blkSize, Exceptions::RuntimeError, "BlockedCoordinatesTransferFactory: Subblocks have different Block sizes. This is not yet supported.");
136 }
137
138 GO indexBase = coarseMap->getFullMap()->getIndexBase();
139 size_t numElements = elementAList.size() / blkSize;
140 Array<GO> elementList(numElements);
141
142 // Amalgamate the map
143 for (LO i = 0; i < Teuchos::as<LO>(numElements); i++)
144 elementList[i] = (elementAList[i*blkSize]-indexBase)/blkSize + indexBase;
145
146 RCP<const Map> coarseCoordMap = MapFactory::Build(coarseMap->getFullMap()->lib(),
147 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), elementList, indexBase, coarseMap->getFullMap()->getComm());
148
149 coarseCoordMapBlocked = rcp(new BlockedMap(coarseCoordMap, subBlockMaps, thyraMode));
150 }
151
152 // Build blocked coordinates vector
153 RCP<dBV> bcoarseCoords = rcp(new dBV(coarseCoordMapBlocked,subBlockCoords));
154
155 // Turn the blocked coordinates vector into an unblocked one
156 RCP<dMV> coarseCoords = bcoarseCoords->Merge();
157 Set<RCP<dMV> >(coarseLevel, "Coordinates", coarseCoords);
158 }
159
160 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
162 subFactories_.push_back(factory);
163 }
164
165
166
167} // namespace MueLu
168
169#endif // MUELU_BLOCKEDCOORDINATESTRANSFER_FACTORY_DEF_HPP
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void DeclareInput(Level &finelevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
void AddFactory(const RCP< const FactoryBase > &factory)
Add (sub) coords factory in the end of list of factories in BlockedCoordinatesTransferFactory.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
Namespace for MueLu classes and methods.
@ Runtime0
One-liner description of what is happening.
@ Runtime1
Description of what is happening (more verbose)