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