C++ AssemblyOptions::acaOptions方法代码示例

std::pair<
shared_ptr<typename HypersingularIntegralOperator<
BasisFunctionType, KernelType, ResultType>::LocalAssembler>,
shared_ptr<typename HypersingularIntegralOperator<
BasisFunctionType, KernelType, ResultType>::LocalAssembler>
>
HypersingularIntegralOperator<BasisFunctionType, KernelType, ResultType>::makeAssemblers(
        const QuadratureStrategy& quadStrategy,
        const AssemblyOptions& options) const
{
    typedef Fiber::RawGridGeometry<CoordinateType> RawGridGeometry;
    typedef std::vector<const Fiber::Shapeset<BasisFunctionType>*> ShapesetPtrVector;

    const bool verbose = (options.verbosityLevel() >= VerbosityLevel::DEFAULT);

    shared_ptr<RawGridGeometry> testRawGeometry, trialRawGeometry;
    shared_ptr<GeometryFactory> testGeometryFactory, trialGeometryFactory;
    shared_ptr<Fiber::OpenClHandler> openClHandler;
    shared_ptr<ShapesetPtrVector> testShapesets, trialShapesets;
    bool cacheSingularIntegrals;

    if (verbose)
        std::cout << "Collecting data for assembler construction..." << std::endl;
       this->collectDataForAssemblerConstruction(options,
                                        testRawGeometry, trialRawGeometry,
                                        testGeometryFactory, trialGeometryFactory,
                                        testShapesets, trialShapesets,
                                        openClHandler, cacheSingularIntegrals);
    if (verbose)
        std::cout << "Data collection finished." << std::endl;

    bool makeSeparateOffDiagonalAssembler =
        options.assemblyMode() == AssemblyOptions::ACA &&
        options.acaOptions().mode == AcaOptions::HYBRID_ASSEMBLY;

    return reallyMakeAssemblers(quadStrategy,
                                testGeometryFactory, trialGeometryFactory,
                                testRawGeometry, trialRawGeometry,
                                testShapesets, trialShapesets, openClHandler,
                                options.parallelizationOptions(),
                                options.verbosityLevel(),
                                cacheSingularIntegrals,
                                makeSeparateOffDiagonalAssembler);
}

//.........这里部分代码省略.........

  // Estimate number of entries in each row

  //    This will be useful when we begin to use MPI
  //    // Get global DOF indices for which this process is responsible
  //    const int testGlobalDofCount = testSpace.globalDofCount();
  //    Epetra_Map rowMap(testGlobalDofCount, 0 /* index-base */, comm);
  //    std::vector<int> myTestGlobalDofs(rowMap.MyGlobalElements(),
  //                                      rowMap.MyGlobalElements() +
  //                                      rowMap.NumMyElements());
  //    const int myTestGlobalDofCount = myTestGlobalDofs.size();

  const int testGlobalDofCount = testSpace.globalDofCount();
  const int trialGlobalDofCount = trialSpace.globalDofCount();
  arma::Col<int> nonzeroEntryCountEstimates(testGlobalDofCount);
  nonzeroEntryCountEstimates.fill(0);

  // Upper estimate for the number of global trial DOFs coupled to a given
  // global test DOF: sum of the local trial DOF counts for each element that
  // contributes to the global test DOF in question
  for (size_t e = 0; e < elementCount; ++e)
    for (size_t testLdof = 0; testLdof < testGdofs[e].size(); ++testLdof) {
      int testGdof = testGdofs[e][testLdof];
      if (testGdof >= 0)
        nonzeroEntryCountEstimates(testGdof) += trialGdofs[e].size();
    }

  Epetra_SerialComm comm; // To be replaced once we begin to use MPI
  Epetra_LocalMap rowMap(testGlobalDofCount, 0 /* index_base */, comm);
  Epetra_LocalMap colMap(trialGlobalDofCount, 0 /* index_base */, comm);
  shared_ptr<Epetra_FECrsMatrix> result =
      boost::make_shared<Epetra_FECrsMatrix>(
          Copy, rowMap, colMap, nonzeroEntryCountEstimates.memptr());

  // TODO: make each process responsible for a subset of elements
  // Find maximum number of local dofs per element
  size_t maxLdofCount = 0;
  for (size_t e = 0; e < elementCount; ++e)
    maxLdofCount =
        std::max(maxLdofCount, testGdofs[e].size() * trialGdofs[e].size());

  // Initialise sparse matrix with zeros at required positions
  arma::Col<double> zeros(maxLdofCount);
  zeros.fill(0.);
  for (size_t e = 0; e < elementCount; ++e)
    result->InsertGlobalValues(testGdofs[e].size(), &testGdofs[e][0],
                               trialGdofs[e].size(), &trialGdofs[e][0],
                               zeros.memptr());
  // Add contributions from individual elements
  for (size_t e = 0; e < elementCount; ++e)
    epetraSumIntoGlobalValues(*result, testGdofs[e], trialGdofs[e],
                              localResult[e]);
  result->GlobalAssemble();

  // If assembly mode is equal to ACA and we have AHMED,
  // construct the block cluster tree. Otherwise leave it uninitialized.
  typedef ClusterConstructionHelper<BasisFunctionType> CCH;
  typedef AhmedDofWrapper<CoordinateType> AhmedDofType;
  typedef ExtendedBemCluster<AhmedDofType> AhmedBemCluster;
  typedef bbxbemblcluster<AhmedDofType, AhmedDofType> AhmedBemBlcluster;

  shared_ptr<AhmedBemBlcluster> blockCluster;
  shared_ptr<IndexPermutation> test_o2pPermutation, test_p2oPermutation;
  shared_ptr<IndexPermutation> trial_o2pPermutation, trial_p2oPermutation;
#ifdef WITH_AHMED
  if (options.assemblyMode() == AssemblyOptions::ACA) {
    const AcaOptions &acaOptions = options.acaOptions();
    bool indexWithGlobalDofs = acaOptions.mode != AcaOptions::HYBRID_ASSEMBLY;

    typedef ClusterConstructionHelper<BasisFunctionType> CCH;
    shared_ptr<AhmedBemCluster> testClusterTree;
    CCH::constructBemCluster(testSpace, indexWithGlobalDofs, acaOptions,
                             testClusterTree, test_o2pPermutation,
                             test_p2oPermutation);
    // TODO: construct a hermitian H-matrix if possible
    shared_ptr<AhmedBemCluster> trialClusterTree;
    CCH::constructBemCluster(trialSpace, indexWithGlobalDofs, acaOptions,
                             trialClusterTree, trial_o2pPermutation,
                             trial_p2oPermutation);
    unsigned int blockCount = 0;
    bool useStrongAdmissibilityCondition = !indexWithGlobalDofs;
    blockCluster.reset(CCH::constructBemBlockCluster(
        acaOptions, false /* hermitian */, *testClusterTree, *trialClusterTree,
        useStrongAdmissibilityCondition, blockCount).release());
  }
#endif

  // Create and return a discrete operator represented by the matrix that
  // has just been calculated
  return std::unique_ptr<DiscreteBoundaryOperator<ResultType>>(
      new DiscreteSparseBoundaryOperator<ResultType>(
          result, this->symmetry(), NO_TRANSPOSE, blockCluster,
          trial_o2pPermutation, test_o2pPermutation));
#else // WITH_TRILINOS
  throw std::runtime_error(
      "ElementaryLocalOperator::assembleWeakFormInSparseMode(): "
      "To enable assembly in sparse mode, recompile BEM++ "
      "with the symbol WITH_TRILINOS defined.");
#endif
}

std::auto_ptr<DiscreteBoundaryOperator<ResultType> >
AcaGlobalAssembler<BasisFunctionType, ResultType>::assembleDetachedWeakForm(
        const Space<BasisFunctionType>& testSpace,
        const Space<BasisFunctionType>& trialSpace,
        const std::vector<LocalAssembler*>& localAssemblers,
        const std::vector<const DiscreteBndOp*>& sparseTermsToAdd,
        const std::vector<ResultType>& denseTermsMultipliers,
        const std::vector<ResultType>& sparseTermsMultipliers,
        const AssemblyOptions& options,
        int symmetry)
{
#ifdef WITH_AHMED
    typedef AhmedDofWrapper<CoordinateType> AhmedDofType;
    typedef ExtendedBemCluster<AhmedDofType> AhmedBemCluster;
    typedef bemblcluster<AhmedDofType, AhmedDofType> AhmedBemBlcluster;
    typedef DiscreteAcaBoundaryOperator<ResultType> DiscreteAcaLinOp;

    const AcaOptions& acaOptions = options.acaOptions();
    const bool indexWithGlobalDofs = acaOptions.globalAssemblyBeforeCompression;
    const bool verbosityAtLeastDefault =
            (options.verbosityLevel() >= VerbosityLevel::DEFAULT);
    const bool verbosityAtLeastHigh =
            (options.verbosityLevel() >= VerbosityLevel::HIGH);

    // Currently we don't support Hermitian ACA operators. This is because we
    // don't have the means to really test them -- we would need complex-valued
    // basis functions for that. (Assembly of such a matrix would be very easy
    // -- just change complex_sym from true to false in the call to apprx_sym()
    // in AcaWeakFormAssemblerLoopBody::operator() -- but operations on
    // symmetric/Hermitian matrices are not always trivial and we do need to be
    // able to test them properly.)
    bool symmetric = symmetry & SYMMETRIC;
    if (symmetry & HERMITIAN && !(symmetry & SYMMETRIC) &&
            verbosityAtLeastDefault)
        std::cout << "Warning: assembly of non-symmetric Hermitian H-matrices "
                     "is not supported yet. A general H-matrix will be assembled"
                  << std::endl;

#ifndef WITH_TRILINOS
    if (!indexWithGlobalDofs)
        throw std::runtime_error("AcaGlobalAssembler::assembleDetachedWeakForm(): "
                                 "ACA assembly with globalAssemblyBeforeCompression "
                                 "set to false requires BEM++ to be linked with "
                                 "Trilinos");
#endif // WITH_TRILINOS

    const size_t testDofCount = indexWithGlobalDofs ?
                testSpace.globalDofCount() : testSpace.flatLocalDofCount();
    const size_t trialDofCount = indexWithGlobalDofs ?
                trialSpace.globalDofCount() : trialSpace.flatLocalDofCount();

    if (symmetric && testDofCount != trialDofCount)
        throw std::invalid_argument("AcaGlobalAssembler::assembleDetachedWeakForm(): "
                                    "you cannot generate a symmetric weak form "
                                    "using test and trial spaces with different "
                                    "numbers of DOFs");

    // o2p: map of original indices to permuted indices
    // p2o: map of permuted indices to original indices
    typedef ClusterConstructionHelper<BasisFunctionType> CCH;
    shared_ptr<AhmedBemCluster> testClusterTree;
    shared_ptr<IndexPermutation> test_o2pPermutation, test_p2oPermutation;
    CCH::constructBemCluster(testSpace, indexWithGlobalDofs, acaOptions,
                             testClusterTree,
                             test_o2pPermutation, test_p2oPermutation);
    shared_ptr<AhmedBemCluster> trialClusterTree;
    shared_ptr<IndexPermutation> trial_o2pPermutation, trial_p2oPermutation;
    if (symmetric || &testSpace == &trialSpace) {
        trialClusterTree = testClusterTree;
        trial_o2pPermutation = test_o2pPermutation;
        trial_p2oPermutation = test_p2oPermutation;
    } else
        CCH::constructBemCluster(trialSpace, indexWithGlobalDofs, acaOptions,
                                 trialClusterTree,
                                 trial_o2pPermutation, trial_p2oPermutation);

//    // Export VTK plots showing the disctribution of leaf cluster ids
//    std::vector<unsigned int> testClusterIds;
//    getClusterIds(*testClusterTree, test_p2oPermutation->permutedIndices(), testClusterIds);
//    testSpace.dumpClusterIds("testClusterIds", testClusterIds,
//                             indexWithGlobalDofs ? GLOBAL_DOFS : FLAT_LOCAL_DOFS);
//    std::vector<unsigned int> trialClusterIds;
//    getClusterIds(*trialClusterTree, trial_p2oPermutation->permutedIndices(), trialClusterIds);
//    trialSpace.dumpClusterIds("trialClusterIds", trialClusterIds,
//                              indexWithGlobalDofs ? GLOBAL_DOFS : FLAT_LOCAL_DOFS);

    if (verbosityAtLeastHigh)
        std::cout << "Test cluster count: " << testClusterTree->getncl()
                  << "\nTrial cluster count: " << trialClusterTree->getncl()
                  << std::endl;

    unsigned int blockCount = 0;
    shared_ptr<AhmedBemBlcluster> bemBlclusterTree(
                CCH::constructBemBlockCluster(acaOptions, symmetric,
                                              *testClusterTree, *trialClusterTree,
                                              blockCount).release());

    if (verbosityAtLeastHigh)
        std::cout << "Mblock count: " << blockCount << std::endl;

//.........这里部分代码省略.........