C++ Basis::evaluate方法代码示例

void NumericalTestTrialIntegrator<BasisFunctionType, ResultType, GeometryFactory>::integrate(
        const std::vector<int>& elementIndices,
        const Basis<BasisFunctionType>& testBasis,
        const Basis<BasisFunctionType>& trialBasis,
        arma::Cube<ResultType>& result) const
{
    const size_t pointCount = m_localQuadPoints.n_cols;
    const size_t elementCount = elementIndices.size();

    if (pointCount == 0 || elementCount == 0)
        return;
    // TODO: in the (pathological) case that pointCount == 0 but
    // elementCount != 0, set elements of result to 0.

    // Evaluate constants
    const int componentCount = m_testTransformations.resultDimension(0);
    const int testDofCount = testBasis.size();
    const int trialDofCount = trialBasis.size();

//    if (m_trialTransformations.codomainDimension() != componentCount)
//        throw std::runtime_error("NumericalTestTrialIntegrator::integrate(): "
//                                 "test and trial functions "
//                                 "must have the same number of components");

    BasisData<BasisFunctionType> testBasisData, trialBasisData;
    GeometricalData<CoordinateType> geomData;

    size_t testBasisDeps = 0, trialBasisDeps = 0;
    size_t geomDeps = INTEGRATION_ELEMENTS;

    m_testTransformations.addDependencies(testBasisDeps, geomDeps);
    m_trialTransformations.addDependencies(trialBasisDeps, geomDeps);

    typedef typename GeometryFactory::Geometry Geometry;
    std::auto_ptr<Geometry> geometry(m_geometryFactory.make());

    CollectionOf3dArrays<BasisFunctionType> testValues, trialValues;

    result.set_size(testDofCount, trialDofCount, elementCount);

    testBasis.evaluate(testBasisDeps, m_localQuadPoints, ALL_DOFS, testBasisData);
    trialBasis.evaluate(trialBasisDeps, m_localQuadPoints, ALL_DOFS, trialBasisData);

    // Iterate over the elements
    for (size_t e = 0; e < elementCount; ++e)
    {
        m_rawGeometry.setupGeometry(elementIndices[e], *geometry);
        geometry->getData(geomDeps, m_localQuadPoints, geomData);
        m_testTransformations.evaluate(testBasisData, geomData, testValues);
        m_trialTransformations.evaluate(trialBasisData, geomData, trialValues);

        for (int trialDof = 0; trialDof < trialDofCount; ++trialDof)
            for (int testDof = 0; testDof < testDofCount; ++testDof)
            {
                ResultType sum = 0.;
                for (size_t point = 0; point < pointCount; ++point)
                    for (int dim = 0; dim < componentCount; ++dim)
                        sum +=  m_quadWeights[point] *
                                geomData.integrationElements(point) *
                                conjugate(testValues[0](dim, testDof, point)) *
                                trialValues[0](dim, trialDof, point);
                result(testDof, trialDof, e) = sum;
            }
    }
}

BOOST_AUTO_TEST_CASE_TEMPLATE(works, BasisFunctionType, basis_function_types)
{
    typedef BasisFunctionType BFT;
    typedef typename Fiber::ScalarTraits<BFT>::ComplexType RT;
    typedef typename Fiber::ScalarTraits<BFT>::RealType CT;
    GridParameters params;
    params.topology = GridParameters::TRIANGULAR;
    shared_ptr<Grid> grid = GridFactory::importGmshGrid(
                params, "meshes/two_disjoint_triangles.msh",
                false /* verbose */);

    PiecewiseLinearContinuousScalarSpace<BFT> pwiseLinears(grid);
    PiecewiseConstantScalarSpace<BFT> pwiseConstants(grid);

    AssemblyOptions assemblyOptions;
    assemblyOptions.setVerbosityLevel(VerbosityLevel::LOW);
    AccuracyOptions accuracyOptions;
    accuracyOptions.doubleRegular.setAbsoluteQuadratureOrder(5);
    accuracyOptions.doubleSingular.setAbsoluteQuadratureOrder(5);
    NumericalQuadratureStrategy<BFT, RT> quadStrategy(accuracyOptions);

    Context<BFT, RT> context(make_shared_from_ref(quadStrategy), assemblyOptions);

    const RT waveNumber(1.23, 0.31);

    BoundaryOperator<BFT, RT> slpOpConstants = Bempp::helmholtz3dSingleLayerBoundaryOperator<BFT>(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseConstants),
                make_shared_from_ref(pwiseConstants),
                make_shared_from_ref(pwiseConstants),
                waveNumber);
    BoundaryOperator<BFT, RT> slpOpLinears = helmholtz3dSingleLayerBoundaryOperator<BFT>(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                waveNumber);
    BoundaryOperator<BFT, RT> hypOp = helmholtz3dHypersingularBoundaryOperator<BFT>(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                waveNumber);

    // Get the matrix repr. of the hypersingular operator
    arma::Mat<RT> hypMat = hypOp.weakForm()->asMatrix();

    // Construct the expected hypersingular operator matrix. For this, we need:

    // * the surface curls of all basis functions (which are constant)
    typedef Fiber::SurfaceCurl3dElementaryFunctor<CT> ElementaryFunctor;
    typedef Fiber::ElementaryBasisTransformationFunctorWrapper<ElementaryFunctor> Functor;
    Functor functor;
    size_t basisDeps = 0, geomDeps = 0;
    functor.addDependencies(basisDeps, geomDeps);

    arma::Mat<CT> points(2, 1);
    points.fill(0.);

    typedef Fiber::PiecewiseLinearContinuousScalarBasis<3, BFT> Basis;
    Basis basis;
    Fiber::BasisData<BFT> basisData;
    basis.evaluate(basisDeps, points, Fiber::ALL_DOFS, basisData);

    Fiber::GeometricalData<CT> geomData[2];
    std::unique_ptr<GridView> view = grid->leafView();
    std::unique_ptr<EntityIterator<0> > it = view->entityIterator<0>();
    it->entity().geometry().getData(geomDeps, points, geomData[0]);
    it->next();
    it->entity().geometry().getData(geomDeps, points, geomData[1]);

    Fiber::DefaultCollectionOfBasisTransformations<Functor> transformations(functor);

    Fiber::CollectionOf3dArrays<BFT> surfaceCurl[2];
    transformations.evaluate(basisData, geomData[0], surfaceCurl[0]);
    transformations.evaluate(basisData, geomData[1], surfaceCurl[1]);

    // * the single-layer potential matrix for constant basis functions
    arma::Mat<RT> slpMatConstants = slpOpConstants.weakForm()->asMatrix();
    // * the single-layer potential matrix for linear basis functions
    arma::Mat<RT> slpMatLinears = slpOpLinears.weakForm()->asMatrix();

    arma::Mat<RT> expectedHypMat(6, 6);
    for (size_t testElement = 0; testElement < 2; ++testElement)
        for (size_t trialElement = 0; trialElement < 2; ++trialElement)
            for (size_t r = 0; r < 3; ++r)
                for (size_t c = 0; c < 3; ++c) {
                    RT curlMultiplier = 0.;
                    for (size_t dim = 0; dim < 3; ++dim)
                        curlMultiplier += surfaceCurl[testElement][0](dim, r, 0) *
                                surfaceCurl[trialElement][0](dim, c, 0);
                    RT valueMultiplier = 0.;
                    for (size_t dim = 0; dim < 3; ++dim)
                        valueMultiplier += geomData[testElement].normals(dim, 0) *
                            geomData[trialElement].normals(dim, 0);
                    expectedHypMat(3 * testElement + r, 3 * trialElement + c) =
                            curlMultiplier * slpMatConstants(testElement, trialElement) -
                            waveNumber * waveNumber * valueMultiplier *
                            slpMatLinears(3 * testElement + r, 3 * trialElement + c);
                }
//.........这里部分代码省略.........

int main(int argc, char** argv) {
    
    const int d  = 3;
    const int dT = 3;
    const int jmax = 10;
    const bool normalization = 0;//choose 1 for Laplacian
    
    const unsigned int testcase=5;
    PeriodicTestProblem<testcase> tper;
    Function<1>* uexact = 0;
    switch(testcase) {
        case 1:
            uexact = new Function1();
            break;
        case 2:
            uexact = new Function2();
            break;
        case 3:
            uexact = new Function3();
            break;
        case 4:
            uexact = new Function4();
            break;
        case 5:
            uexact = new Hat();    
        default:
            break;
    }
    
    
    
#ifdef PERIODIC_CDFBASIS
    typedef CDFBasis<d,dT> RBasis;
    typedef PeriodicBasis<RBasis> Basis;
    typedef Basis::Index Index;
    Basis basis;
    basis.set_jmax(jmax);
    
    typedef PeriodicIntervalGramian<RBasis> Problem;
    Problem G(tper, basis);
    const int j0= G.basis().j0();
    
#if 1
    //Plot of wavelet with coefficient mu
    Index mu(3,0,1);
    SampledMapping<1> sm1(basis.evaluate(mu, 8, normalization));
    std::ofstream u_stream1("plot_periodicwavelet.m");
    sm1.matlab_output(u_stream1);
    u_stream1 << "figure;\nplot(x,y);"
            << "title('periodic wavelet/generator with index: " << mu << "');" << endl;
    u_stream1.close();
#endif
    
    set<Index> Lambda;
    Vector<double> x;
    
  for (int j = j0; j <= jmax; j++) {
        Lambda.clear();
        
    cout << "  j=" << j << ":" << endl;
    //Implementation of the index set
    
    for (Index lambda = G.basis().first_generator(j0);; ++lambda) {
        Lambda.insert(lambda);
        if (lambda == G.basis().last_wavelet(j)) break;    
    }
    SparseMatrix<double> A;
    cout << "- set up stiffness matrix..." << endl;
    setup_stiffness_matrix(G, Lambda, A);
//    
    
    
    cout << "- set up right-hand side..." << endl;
    Vector<double> b;
    setup_righthand_side(G, Lambda, b);
//    cout << "- right hand side: " << b << endl << endl;
    
    
    x.resize(Lambda.size()); x = 0;
    Vector<double> residual(Lambda.size()); 
    unsigned int iterations;
    
    CG(A, b, x, 1e-15, 250, iterations);
    cout << "  Galerkin system solved with residual (infinity) norm ";
    A.apply(x, residual);
    residual -= b;
    cout << linfty_norm(residual)
	 << " (" << iterations << " iterations needed)" << endl;
    
    // evaluate approximate solution on a grid
    const unsigned int N = 100;
    const double h = 1./N;
    Vector<double> u_values(N+1);
    for (unsigned int i = 0; i <= N; i++) {
      const double point = i*h;
      int id = 0;
      for (set<Index>::const_iterator it(Lambda.begin()); it != Lambda.end(); ++it, ++id) 
	u_values[i] += x[id] * basis.evaluate(0, *it, point, normalization);
    }
//     cout << "  point values of Galerkin solution: " << u_values << endl;
//.........这里部分代码省略.........