Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_AdaptivityUtils.cpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
47
48#include "Epetra_IntVector.h"
49#include "Epetra_Import.h"
50
52 const Epetra_Comm & Comm,
53 const std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > & per_dof_row_basis,
54 std::vector<int> & myRowGidOffsets)
55{
56 // build row offsets
57 int totalSz = 0;
58 std::vector<int> myRowOffsets(per_dof_row_basis.size());
59 for(std::size_t i=0;i<per_dof_row_basis.size();i++) {
60 myRowOffsets[i] = totalSz;
61 totalSz += per_dof_row_basis[i]->size();
62 }
63
64 // build row map
65 Teuchos::RCP<Epetra_Map> rowMap = Teuchos::rcp(new Epetra_Map(-1,totalSz,0,Comm));
66
67 myRowGidOffsets.resize(per_dof_row_basis.size());
68 for(std::size_t i=0;i<myRowGidOffsets.size();i++)
69 myRowGidOffsets[i] = rowMap->GID(myRowOffsets[i]);
70
71 return rowMap;
72}
73
75 const Epetra_Comm & Comm,
76 const std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > & per_dof_row_basis)
77{
78 // build row offsets
79 int totalSz = 0;
80 for(std::size_t i=0;i<per_dof_row_basis.size();i++)
81 totalSz += per_dof_row_basis[i]->size();
82
83 // build row map
84 return Teuchos::rcp(new Epetra_Map(-1,totalSz,0,Comm));
85}
86
88 const Epetra_CrsGraph & determGraph,
89 const std::vector<int> & myRowGidOffsets,
90 std::vector<int> & myColGidOffsets)
91{
92 // build global distributed offsets index
93 Epetra_IntVector adaptRowGids(determGraph.RowMap()); // uniquely defined row GIDs that will
94 // be used to determine col GIDs
95 for(std::size_t i=0;i<myRowGidOffsets.size();i++)
96 adaptRowGids[i] = myRowGidOffsets[i];
97
98 // perform global communication to determine offsets
99 Epetra_Import importer(determGraph.ColMap(),determGraph.RowMap());
100 Epetra_IntVector adaptColGids(determGraph.ColMap());
101 adaptColGids.Import(adaptRowGids,importer,Insert);
102
103 // copy offsets into std::vector
104 myColGidOffsets.resize(adaptColGids.MyLength());
105 for(std::size_t i=0;i<myColGidOffsets.size();i++)
106 myColGidOffsets[i] = adaptColGids[i];
107}
108
110 const Epetra_CrsGraph & determGraph,
111 const Teuchos::RCP<const Stokhos::ProductBasis<int,double> > & masterBasis,
112 const std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > & per_dof_row_basis,
113 std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > & per_dof_col_basis)
114{
115 using Teuchos::RCP;
116
117 const Epetra_BlockMap & rowMap = determGraph.RowMap();
118 const Epetra_BlockMap & colMap = determGraph.ColMap();
119
120 // this is block size
121 int numStochDim = masterBasis->dimension();
122
123 // build blocked maps
124 Epetra_BlockMap blkRowMap(-1,rowMap.NumMyElements(),rowMap.MyGlobalElements(),numStochDim,0,rowMap.Comm());
125 Epetra_BlockMap blkColMap(-1,colMap.NumMyElements(),colMap.MyGlobalElements(),numStochDim,0,colMap.Comm());
126 // construct int vectors to determine order
127 Epetra_IntVector stochRowOrders(blkRowMap),
128 stochColOrders(blkColMap); // to build built by Import
129
130 // loop over row degrees of freedom building Row Orders vector
131 for(std::size_t dof=0;dof<per_dof_row_basis.size();dof++) {
132 RCP<const Stokhos::ProductBasis<int,double> > rowBasis = per_dof_row_basis[dof];
133 TEUCHOS_TEST_FOR_EXCEPTION(rowBasis->dimension()!=masterBasis->dimension(),std::invalid_argument,
134 "Stokhos::adapt_utils::buildColBasisFunctions: Row basis must match dimension of master basis!");
135
136 Teuchos::Array<RCP<const OneDOrthogPolyBasis<int,double> > > onedBasis
137 = rowBasis->getCoordinateBases();
138
139 TEUCHOS_TEST_FOR_EXCEPTION(onedBasis.size()!=numStochDim,std::logic_error,
140 "Stokhos::adapt_utils::buildColBasisFunctions: Wrong number of dimensions from row basis!");
141
142 // fill stochastic orders vector
143 for(int i=0;i<numStochDim;i++)
144 stochRowOrders[i+dof*numStochDim] = onedBasis[i]->order();
145 }
146
147 // do communication to determine row maps
148 Epetra_Import importer(blkColMap,blkRowMap);
149 stochColOrders.Import(stochRowOrders,importer,Insert);
150
151 Teuchos::Array<RCP<const OneDOrthogPolyBasis<int,double> > > oneDMasterBasis
152 = masterBasis->getCoordinateBases();
153
154 // now construct column basis functions
155 std::vector<int> polyOrder(numStochDim,0);
156 per_dof_col_basis.resize(blkColMap.NumMyElements());
157 for(std::size_t col=0;col<per_dof_col_basis.size();col++) {
158 int colOffset = numStochDim*col;
159
160 // extract polynomial order for this column
161 for(int o=0;o<numStochDim;o++)
162 polyOrder[o] = stochColOrders[colOffset+o];
163
164 Teuchos::Array<Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > newBasisArray(numStochDim);
165 for(Teuchos::Ordinal dim=0;dim<oneDMasterBasis.size();dim++) {
166 RCP<const OneDOrthogPolyBasis<int,double> > oneDBasis = oneDMasterBasis[dim];
167
168 newBasisArray[dim] = oneDBasis->cloneWithOrder(polyOrder[dim]);
169 }
170
171 per_dof_col_basis[col] = Teuchos::rcp(
173 }
174}
175
176Teuchos::RCP<Epetra_CrsGraph> Stokhos::adapt_utils::buildAdaptedGraph(
177 const Epetra_CrsGraph & determGraph,
178 const Teuchos::RCP<const Stokhos::ProductBasis<int,double> > & masterBasis,
179 const std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > & per_dof_row_basis,
180 bool onlyUseLinear,
181 int kExpOrder)
182{
183 std::vector<int> myRowGidOffsets, myColGidOffsets;
184
185 return buildAdaptedGraph(determGraph,masterBasis,per_dof_row_basis,myRowGidOffsets,myColGidOffsets,onlyUseLinear,kExpOrder);
186}
187
188Teuchos::RCP<Epetra_CrsGraph> Stokhos::adapt_utils::buildAdaptedGraph(
189 const Epetra_CrsGraph & determGraph,
190 const Teuchos::RCP<const Stokhos::ProductBasis<int,double> > & masterBasis,
191 const std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > & per_dof_row_basis,
192 std::vector<int> & myRowGidOffsets,std::vector<int> & myColGidOffsets,
193 bool onlyUseLinear,
194 int kExpOrder)
195{
196 TEUCHOS_TEST_FOR_EXCEPTION(int(per_dof_row_basis.size())!=determGraph.NumMyRows(),std::logic_error,
197 "Stokhos::adapt_utils::buildAdaptedGraph: per_dof_row_basis.size()!=determGraph.NumMyRows()");
198
199 myRowGidOffsets.clear();
200 myColGidOffsets.clear();
201
202 std::vector<Teuchos::RCP<const Stokhos::ProductBasis<int,double> > > per_dof_col_basis;
203 Teuchos::RCP<Epetra_Map> rowMap;
204
205 // build basis functions associated with the columns
206 buildColBasisFunctions(determGraph,masterBasis,per_dof_row_basis,per_dof_col_basis);
207
208 // build row map, and row and column offsets.
209 rowMap = buildAdaptedRowMapAndOffsets(determGraph.Comm(),per_dof_row_basis,myRowGidOffsets);
210 buildAdaptedColOffsets(determGraph,myRowGidOffsets,myColGidOffsets);
211
212 Teuchos::RCP<Epetra_CrsGraph> graph = Teuchos::rcp(new Epetra_CrsGraph(Copy,*rowMap,0));
213
214 Teuchos::RCP<const Stokhos::Sparse3Tensor<int,double> > Cijk;
215 if(kExpOrder<0)
216 Cijk = masterBasis->computeTripleProductTensor();
217 else
218 Cijk = masterBasis->computeLinearTripleProductTensor();
219
220 // iterate over nonzero structure of graph
221 int maxNNZ = determGraph.MaxNumNonzeros();
222 std::vector<int> determGraphCols(maxNNZ);
223 std::vector<int> graphCols;
224 for(int lRID=0;lRID<determGraph.NumMyRows();lRID++) {
225 int gRID = determGraph.GRID(lRID);
226 int numIndices = -1;
227 int rowOffsetIndex = myRowGidOffsets[lRID];
228
229 // grab row nonzero structure
230 determGraph.ExtractGlobalRowCopy(gRID,maxNNZ,numIndices,&determGraphCols[0]);
231
232 for(int i=0;i<numIndices;i++) {
233 int gCID = determGraphCols[i];
234 int lCID = determGraph.LCID(gCID);
235 int colOffsetIndex = myColGidOffsets[lCID];
236
237 Stokhos::BasisInteractionGraph interactGraph(*masterBasis,*Cijk,*per_dof_row_basis[lRID],
238 *per_dof_col_basis[lCID],
239 onlyUseLinear,kExpOrder);
240 for(std::size_t basisRow=0;basisRow<interactGraph.rowCount();basisRow++) {
241 const std::vector<std::size_t> & basisCols = interactGraph.activeIndices(basisRow);
242 graphCols.resize(basisCols.size());
243 for(std::size_t g=0;g<basisCols.size();g++)
244 graphCols[g] = basisCols[g] + colOffsetIndex;
245
246 graph->InsertGlobalIndices(basisRow+rowOffsetIndex,graphCols.size(),&graphCols[0]);
247 }
248 }
249 }
250
251 graph->FillComplete();
252
253 return graph;
254}
Insert
Copy
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
Teuchos::RCP< Epetra_Map > buildAdaptedRowMapAndOffsets(const Epetra_Comm &Comm, const std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &per_dof_row_basis, std::vector< int > &myRowGidOffsets)
void buildColBasisFunctions(const Epetra_CrsGraph &determGraph, const Teuchos::RCP< const Stokhos::ProductBasis< int, double > > &masterBasis, const std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &per_dof_row_basis, std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &per_dof_col_basis)
Teuchos::RCP< Epetra_Map > buildAdaptedRowMap(const Epetra_Comm &Comm, const std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &per_dof_row_basis)
Teuchos::RCP< Epetra_CrsGraph > buildAdaptedGraph(const Epetra_CrsGraph &determGraph, const Teuchos::RCP< const Stokhos::ProductBasis< int, double > > &masterBasis, const std::vector< Teuchos::RCP< const Stokhos::ProductBasis< int, double > > > &per_dof_row_basis, bool onlyUseLinear=false, int kExpOrder=-1)
void buildAdaptedColOffsets(const Epetra_CrsGraph &determGraph, const std::vector< int > &myRowGidOffsets, std::vector< int > &myColGidOffsets)