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 "SundanceVectorFillingAssemblyKernel.hpp"
00034 #include "Teuchos_Time.hpp"
00035 #include "Teuchos_TimeMonitor.hpp"
00036
00037
00038 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00039 #include "PlayaVectorImpl.hpp"
00040 #endif
00041
00042 using namespace Sundance;
00043 using namespace Sundance;
00044 using namespace Sundance;
00045 using namespace Teuchos;
00046 using namespace Playa;
00047 using std::setw;
00048 using std::endl;
00049
00050 static Time& vecInsertTimer()
00051 {
00052 static RCP<Time> rtn
00053 = TimeMonitor::getNewTimer("vector insertion");
00054 return *rtn;
00055 }
00056
00057 VectorFillingAssemblyKernel::VectorFillingAssemblyKernel(
00058 const Array<RCP<DOFMapBase> >& dofMap,
00059 const Array<RCP<Array<int> > >& isBCIndex,
00060 const Array<int>& lowestLocalIndex,
00061 Array<Vector<double> >& b,
00062 bool partitionBCs,
00063 int verbosity
00064 )
00065 : AssemblyKernelBase(verbosity),
00066 b_(b),
00067 vec_(b.size()),
00068 mapBundle_(dofMap, isBCIndex, lowestLocalIndex, partitionBCs, verbosity)
00069 {
00070 Tabs tab0;
00071
00072 SUNDANCE_MSG1(verb(), tab0 << "VectorFillingAssemblyKernel ctor");
00073
00074 int numBlocks = dofMap.size();
00075
00076 for (int i=0; i<b_.size(); i++)
00077 {
00078 vec_[i].resize(numBlocks);
00079 for (int block=0; block<numBlocks; block++)
00080 {
00081 Tabs tab1;
00082 SUNDANCE_MSG1(verb(), tab1 << "getting vector for block b="
00083 << block << " of " << numBlocks);
00084 Vector<double> vecBlock = b[i].getBlock(block);
00085 vec_[i][block] = rcp_dynamic_cast<LoadableVector<double> >(vecBlock.ptr());
00086
00087 TEUCHOS_TEST_FOR_EXCEPTION(vec_[i][block].get()==0, std::runtime_error,
00088 "vector block " << block << " is not loadable");
00089 vecBlock.zero();
00090 }
00091 }
00092 SUNDANCE_MSG1(verb(), tab0 << "done VectorFillingAssemblyKernel ctor");
00093 }
00094
00095
00096 void VectorFillingAssemblyKernel::buildLocalDOFMaps(
00097 const RCP<StdFwkEvalMediator>& mediator,
00098 IntegrationCellSpecifier intCellSpec,
00099 const Array<Set<int> >& requiredFuncs)
00100 {
00101 mapBundle_.buildLocalDOFMaps(mediator, intCellSpec, requiredFuncs,
00102 verb());
00103 }
00104
00105
00106 void VectorFillingAssemblyKernel::insertLocalVectorBatch(
00107 bool isBCRqc,
00108 bool useCofacetCells,
00109 const Array<int>& funcID,
00110 const Array<int>& funcBlock,
00111 const Array<int>& mvIndices,
00112 const Array<double>& localValues) const
00113 {
00114 TimeMonitor timer(vecInsertTimer());
00115 Tabs tab0;
00116
00117 SUNDANCE_MSG1(verb(), tab0 << "inserting local vector batch");
00118 SUNDANCE_MSG4(verb(), tab0 << "vector values are " << localValues);
00119
00120 const MapBundle& mb = mapBundle_;
00121 int nCells = mb.nCells();
00122
00123 for (int i=0; i<funcID.size(); i++)
00124 {
00125 Tabs tab1;
00126 SUNDANCE_MSG2(verb(), tab1 << "function ID = "<< funcID[i]
00127 << " of " << funcID.size());
00128 SUNDANCE_MSG2(verb(), tab1 << "is BC eqn = " << isBCRqc);
00129 SUNDANCE_MSG2(verb(), tab1 << "num cells = " << nCells);
00130 SUNDANCE_MSG2(verb(), tab1 << "using cofacet cells = " << useCofacetCells);
00131 SUNDANCE_MSG2(verb(), tab1 << "multivector index = "
00132 << mvIndices[i]);
00133
00134
00135
00136 int block = funcBlock[i];
00137
00138 const RCP<DOFMapBase>& dofMap = mb.dofMap(block);
00139 int lowestLocalRow = mb.lowestLocalIndex(block);
00140
00141 int chunk = mb.mapStruct(block, useCofacetCells)->chunkForFuncID(funcID[i]);
00142 SUNDANCE_MSG2(verb(), tab1 << "chunk = " << chunk);
00143
00144 int funcIndex = mb.mapStruct(block, useCofacetCells)->indexForFuncID(funcID[i]);
00145 SUNDANCE_MSG2(verb(), tab1 << "func offset into local DOF map = "
00146 << funcIndex);
00147
00148 const Array<int>& dofs = mb.localDOFs(block, useCofacetCells, chunk);
00149 SUNDANCE_MSG4(verb(), tab1 << "local dofs = " << dofs);
00150
00151 int nFuncs = mb.mapStruct(block, useCofacetCells)->numFuncs(chunk);
00152 SUNDANCE_MSG2(verb(), tab1 << "num funcs in chunk = " << nFuncs);
00153
00154 int nNodes = mb.nNodesInChunk(block, useCofacetCells, chunk);
00155 SUNDANCE_MSG2(verb(), tab1 << "num nodes in chunk = " << nNodes);
00156
00157 const Array<int>& isBCIndex = *(mb.isBCIndex(block));
00158
00159
00160 int r=0;
00161 RCP<Playa::LoadableVector<double> > vecBlock
00162 = vec_[mvIndices[i]][block];
00163
00164 FancyOStream& os = Out::os();
00165
00166 for (int c=0; c<nCells; c++)
00167 {
00168 Tabs tab2;
00169 SUNDANCE_MSG2(verb(), tab2 << "cell = " << c << " of " << nCells);
00170 for (int n=0; n<nNodes; n++, r++)
00171 {
00172 Tabs tab3;
00173 int rowIndex = dofs[(c*nFuncs + funcIndex)*nNodes + n];
00174 int localRowIndex = rowIndex - lowestLocalRow;
00175 if (verb() >= 2) os << tab3 << "n=" << setw(4) << n
00176 << " G=" << setw(8) << rowIndex
00177 << " L=" << setw(8) << localRowIndex;
00178 if (!(dofMap->isLocalDOF(rowIndex))
00179 || isBCRqc!=isBCIndex[localRowIndex])
00180 {
00181 if (verb() >= 2)
00182 {
00183 if (!dofMap->isLocalDOF(rowIndex))
00184 {
00185 os << " --- skipping (is non-local) ---" << std::endl;
00186 }
00187 else if (!isBCRqc && isBCIndex[localRowIndex])
00188 {
00189 os << " --- skipping (is BC row) ---" << std::endl;
00190 }
00191 else
00192 {
00193 os << " --- skipping (is non-BC row) ---" << std::endl;
00194 }
00195 }
00196 }
00197 else
00198 {
00199 if (verb() >= 2) os << setw(15) << localValues[r] << std::endl;
00200 vecBlock->addToElement(rowIndex, localValues[r]);
00201 }
00202 }
00203 }
00204 }
00205 SUNDANCE_MSG1(verb(), tab0 << "...done vector insertion");
00206 }
00207
00208