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
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "SundanceOut.hpp"
00043 #include "PlayaTabs.hpp"
00044 #include "SundanceMatrixVectorAssemblyKernel.hpp"
00045 #include "PlayaSimpleZeroOpDecl.hpp"
00046
00047
00048
00049
00050
00051
00052 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00053 #include "PlayaSimpleZeroOpImpl.hpp"
00054 #include "PlayaVectorImpl.hpp"
00055 #endif
00056
00057 using namespace Sundance;
00058 using namespace Sundance;
00059 using namespace Sundance;
00060 using namespace Teuchos;
00061 using namespace Playa;
00062 using std::setw;
00063 using std::setprecision;
00064 using std::ios_base;
00065 using std::endl;
00066
00067
00068
00069 void MatrixVectorAssemblyKernel::init(
00070 const Array<RCP<DOFMapBase> >& rowMap,
00071 const Array<RCP<DOFMapBase> >& colMap,
00072 LinearOperator<double> A,
00073 bool partitionBCs)
00074 {
00075 Tabs tab;
00076 SUNDANCE_MSG2(verb(), tab << "begin MVAssemblyKernel::init()");
00077
00078 int numRowBlocks = rowMap.size();
00079 int numColBlocks = colMap.size();
00080
00081 for (int br=0; br<numRowBlocks; br++)
00082 {
00083 Vector<double> vecBlock;
00084
00085 mat_[br].resize(numColBlocks);
00086 for (int bc=0; bc<numColBlocks; bc++)
00087 {
00088 Tabs tab1;
00089 LinearOperator<double> matBlock;
00090 if (partitionBCs && numRowBlocks==1 && numColBlocks==1)
00091 {
00092 matBlock = A;
00093 }
00094 else
00095 {
00096 matBlock = A.getBlock(br, bc);
00097 }
00098 if (matBlock.ptr().get() == 0) continue;
00099 const SimpleZeroOp<double>* zp
00100 = dynamic_cast<const SimpleZeroOp<double>*>(matBlock.ptr().get());
00101 if (zp)
00102 {
00103 SUNDANCE_MSG3(verb(), tab1 << "block(br=" << br << ", "
00104 << bc << ") is zero");
00105 TEUCHOS_TEST_FOR_EXCEPTION(numRowBlocks==1 && numColBlocks==1 && zp,
00106 std::runtime_error, "no nonzero block in target matrix");
00107 continue;
00108 }
00109 mat_[br][bc]
00110 = dynamic_cast<Playa::LoadableMatrix<double>* >(matBlock.ptr().get());
00111 TEUCHOS_TEST_FOR_EXCEPTION(mat_[br][bc]==0, std::runtime_error,
00112 "matrix block (" << br << ", " << bc
00113 << ") is not loadable in Assembler::assemble()");
00114 mat_[br][bc]->zero();
00115 }
00116 }
00117 SUNDANCE_MSG2(verb(), tab << "end MVAssemblyKernel::init()");
00118 }
00119
00120
00121
00122 void MatrixVectorAssemblyKernel::fill(
00123 bool isBC,
00124 const IntegralGroup& group,
00125 const RCP<Array<double> >& localValues)
00126 {
00127 Tabs tab0;
00128 SUNDANCE_MSG1(verb(), tab0 << "in MatrixVectorAssemblyKernel::fill()");
00129 SUNDANCE_MSG1(verb(), tab0 << "filling for integral group " << group.derivs());
00130
00131 bool useCofacets = group.usesMaximalCofacets();
00132
00133 if (group.isOneForm())
00134 {
00135 Tabs tab1;
00136 SUNDANCE_MSG2(verb(), tab1 << "group is one form");
00137 insertLocalVectorBatch(isBC, useCofacets,
00138 group.testID(), group.testBlock(), group.mvIndices(),
00139 *localValues);
00140 }
00141 else if (group.isTwoForm())
00142 {
00143 Tabs tab1;
00144 SUNDANCE_MSG2(verb(), tab1 << "group is two form");
00145 insertLocalMatrixBatch(isBC, useCofacets,
00146 group.testID(), group.testBlock(),
00147 group.unkID(), group.unkBlock(),
00148 *localValues);
00149 }
00150 else
00151 {
00152 Tabs tab1;
00153 SUNDANCE_MSG2(verb(), tab1 << "is zero form -- nothing to do here");
00154 }
00155
00156 SUNDANCE_MSG1(verb(), tab0 << "done MatrixVectorAssemblyKernel::fill()");
00157 }
00158
00159
00160 void MatrixVectorAssemblyKernel::prepareForWorkSet(
00161 const Array<Set<int> >& requiredTests,
00162 const Array<Set<int> >& requiredUnks,
00163 RCP<StdFwkEvalMediator> mediator)
00164 {
00165 Tabs tab0;
00166 SUNDANCE_MSG1(verb(), tab0
00167 << "in MatrixVectorAssemblyKernel::prepareForWorkSet()");
00168
00169 IntegrationCellSpecifier intCellSpec = mediator->integrationCellSpec();
00170 SUNDANCE_MSG2(verb(), tab0
00171 << "integration cell specifier is " << intCellSpec);
00172
00173 SUNDANCE_MSG2(verb(), tab0 << "building row DOF maps");
00174 buildLocalDOFMaps(mediator, intCellSpec, requiredTests);
00175
00176 SUNDANCE_MSG2(verb(), tab0 << "building column DOF maps");
00177 cmb_.buildLocalDOFMaps(mediator, intCellSpec, requiredUnks, verb());
00178
00179 SUNDANCE_MSG1(verb(), tab0 << "done MatrixVectorAssemblyKernel::prepareForWorkSet()");
00180 }
00181
00182
00183 void MatrixVectorAssemblyKernel::insertLocalMatrixBatch(
00184 bool isBCRqc,
00185 bool useCofacetCells,
00186 const Array<int>& testID,
00187 const Array<int>& testBlock,
00188 const Array<int>& unkID,
00189 const Array<int>& unkBlock,
00190 const Array<double>& localValues) const
00191 {
00192 Tabs tab;
00193
00194 SUNDANCE_MSG1(verb(), tab << "inserting local matrices");
00195
00196 static Array<int> skipRow;
00197 static Array<int> rows;
00198 static Array<int> cols;
00199
00200 int nCells = rmb().nCells();
00201
00202 for (int t=0; t<testID.size(); t++)
00203 {
00204 Tabs tab1;
00205 int br = testBlock[t];
00206
00207 SUNDANCE_MSG3(verb(), tab1 << "block row = " << br
00208 << tab1 << "function ID = "<< testID[t]
00209 << " of " << testID.size() << std::endl
00210 << tab1 << "is BC eqn = " << isBCRqc << std::endl
00211 << tab1 << "num cells = " << nCells << std::endl
00212 << tab1 << "using cofacet cells = " << useCofacetCells);
00213 SUNDANCE_MSG3(verb(), tab1 << "local values=" << localValues);
00214
00215 const RCP<DOFMapBase>& rowMap = rmb().dofMap(br);
00216 int lowestLocalRow = rmb().lowestLocalIndex(br);
00217 int highestRowIndex = lowestLocalRow + rowMap->numLocalDOFs();
00218 int testChunk = rmb().mapStruct(br,
00219 useCofacetCells)->chunkForFuncID(testID[t]);
00220 int testFuncIndex = rmb().mapStruct(br,
00221 useCofacetCells)->indexForFuncID(testID[t]);
00222 const Array<int>& testDofs = rmb().localDOFs(br, useCofacetCells, testChunk);
00223 int nTestFuncs = rmb().mapStruct(br, useCofacetCells)->numFuncs(testChunk);
00224 int nTestNodes = rmb().nNodesInChunk(br, useCofacetCells, testChunk);
00225
00226 SUNDANCE_MSG3(verb(), tab1 << "lowestLocalRow = " << lowestLocalRow << std::endl
00227 << tab1 << "test chunk = " << testChunk << std::endl
00228 << tab1 << "func offset into local DOF map = "
00229 << testFuncIndex << std::endl
00230 << tab1 << "local test dofs = " << testDofs << std::endl
00231 << tab1 << "num test funcs in chunk = " << nTestFuncs << std::endl
00232 << tab1 << "num test nodes in chunk = " << nTestNodes);
00233
00234 int numRows = nCells * nTestNodes;
00235 SUNDANCE_MSG3(verb(), tab1 << "numRows=" << numRows);
00236 const Array<int>& isBCRow = *(rmb().isBCIndex(br));
00237 rows.resize(numRows);
00238 skipRow.resize(numRows);
00239 int r=0;
00240 for (int c=0; c<nCells; c++)
00241 {
00242 for (int n=0; n<nTestNodes; n++, r++)
00243 {
00244 int row = testDofs[(c*nTestFuncs + testFuncIndex)*nTestNodes + n];
00245 rows[r] = row;
00246 int localRow = rows[r]-lowestLocalRow;
00247 skipRow[r] = row < lowestLocalRow || row >= highestRowIndex
00248 || (isBCRqc && !isBCRow[localRow])
00249 || (!isBCRqc && isBCRow[localRow]);
00250 }
00251 }
00252
00253 SUNDANCE_MSG3(verb(), tab1 << "rows=" << rows);
00254 SUNDANCE_MSG3(verb(), tab1 << "skipRow=" << skipRow);
00255 for (int u=0; u<unkID.size(); u++)
00256 {
00257 Tabs tab2;
00258 int bc = unkBlock[u];
00259
00260 int lowestLocalCol = cmb().lowestLocalIndex(bc);
00261 int unkChunk = cmb().mapStruct(bc,
00262 useCofacetCells)->chunkForFuncID(unkID[u]);
00263 int unkFuncIndex = cmb().mapStruct(bc,
00264 useCofacetCells)->indexForFuncID(unkID[u]);
00265 const Array<int>& unkDofs = cmb().localDOFs(bc, useCofacetCells, unkChunk);
00266 int nUnkFuncs = cmb().mapStruct(bc, useCofacetCells)->numFuncs(unkChunk);
00267 int nUnkNodes = cmb().nNodesInChunk(bc, useCofacetCells, unkChunk);
00268
00269
00270 SUNDANCE_MSG3(verb(), tab2 << "lowestLocalCol = "
00271 << lowestLocalCol << std::endl
00272 << tab2 << "block col = " << bc << std::endl
00273 << tab2 << "unk ID = "<< unkID[t]
00274 << " of " << unkID.size() << std::endl
00275 << tab2 << "unk chunk = " << unkChunk << std::endl
00276 << tab2 << "func offset into local DOF map = "
00277 << unkFuncIndex << std::endl
00278 << tab2 << "local unk dofs = " << unkDofs << std::endl
00279 << tab2 << "num unk funcs in chunk = " << nUnkFuncs << std::endl
00280 << tab2 << "num unk nodes in chunk = " << nUnkNodes);
00281
00282 cols.resize(nCells*nUnkNodes);
00283
00284 int j=0;
00285 for (int c=0; c<nCells; c++)
00286 {
00287 for (int n=0; n<nUnkNodes; n++, j++)
00288 {
00289 cols[j] = unkDofs[(c*nUnkFuncs + unkFuncIndex)*nUnkNodes + n];
00290 }
00291 }
00292
00293 if (verb() >= 3)
00294 {
00295 writeLSMs(br, bc, useCofacetCells,
00296 nTestNodes, nTestFuncs, testFuncIndex, testDofs,
00297 nUnkNodes, nUnkFuncs, unkFuncIndex, unkDofs, localValues);
00298 }
00299
00300 SUNDANCE_MSG2(verb(), tab2 << "calling addToElementBatch()");
00301 TEUCHOS_TEST_FOR_EXCEPT(mat_[br][bc]==0);
00302 mat_[br][bc]->addToElementBatch(numRows,
00303 nTestNodes,
00304 &(rows[0]),
00305 nUnkNodes,
00306 &(cols[0]),
00307 &(localValues[0]),
00308 &(skipRow[0]));
00309 SUNDANCE_MSG2(verb(), tab2 << "done calling addToElementBatch()");
00310 }
00311 }
00312
00313 SUNDANCE_MSG1(verb(), tab << "done inserting local matrices");
00314 }
00315
00316
00317 void MatrixVectorAssemblyKernel::writeLSMs(
00318 int blockRow,
00319 int blockCol,
00320 bool useCofacetCells,
00321 int nTestNodes,
00322 int nTestFuncs,
00323 int testFuncIndex,
00324 const Array<int>& testDofs,
00325 int nUnkNodes,
00326 int nUnkFuncs,
00327 int unkFuncIndex,
00328 const Array<int>& unkDofs,
00329 const Array<double>& localValues) const
00330 {
00331 FancyOStream& os = Out::os();
00332
00333 int nCells = rmb().nCells();
00334
00335 RCP<const Array<int> > workSet = rmb().workSet(blockRow, useCofacetCells);
00336
00337 int lr = 0;
00338
00339 ios_base::fmtflags oldFlags = os.flags();
00340 os.setf(ios_base::right);
00341 os.setf(ios_base::showpoint);
00342
00343 for (int c=0; c<nCells; c++)
00344 {
00345 Tabs tab3;
00346 os << tab3 << std::endl
00347 << tab3 << "c="<< c << ", cellLID=" << (*workSet)[c] << std::endl;
00348
00349 os << tab3 << "num values per cell = " << localValues.size()/nCells
00350 << ", num test nodes=" << nTestNodes << ", num unk nodes="
00351 << nUnkNodes << std::endl;
00352 Array<int> lsmCols(nUnkNodes);
00353 os << tab3 << setw(17);
00354
00355 os << "|";
00356 for (int n=0; n<nUnkNodes; n++)
00357 {
00358 int colDof = unkDofs[(c*nUnkFuncs + unkFuncIndex)*nUnkNodes + n];
00359 lsmCols[n] = colDof;
00360 os << setw(12) << colDof;
00361 }
00362 os << std::endl << tab3 << "---------------------------------------------------------------" << std::endl;
00363
00364
00365 for (int m=0; m<nTestNodes; m++, lr++)
00366 {
00367 int globalRow
00368 = testDofs[(c*nTestFuncs + testFuncIndex)*nTestNodes+m];
00369 os << tab3 << setw(6) << m << setw(10) << globalRow << "|";
00370
00371 for (int n=0; n<nUnkNodes; n++)
00372 {
00373 double Amn = localValues[lr*nUnkNodes + n];
00374 os << setw(12) << setprecision(5) << Amn;
00375 }
00376 os << std::endl;
00377 }
00378 }
00379 os.flags(oldFlags);
00380 }
