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