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

BOOST_AUTO_TEST_CASE_TEMPLATE(symmetric_matches_nonsymmetric_in_aca_mode,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename Fiber::ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    if (boost::is_same<RT, std::complex<float> >::value) {
        // The AHMED support for single-precision complex symmetric matrices
        // is broken
        BOOST_CHECK(true);
        return;
    }

    GridParameters params;
    params.topology = GridParameters::TRIANGULAR;
    shared_ptr<Grid> grid = GridFactory::importGmshGrid(
                                params, "../../examples/meshes/sphere-h-0.4.msh",
                                false /* verbose */);

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

    AssemblyOptions assemblyOptions;
    assemblyOptions.setVerbosityLevel(VerbosityLevel::LOW);
    AcaOptions acaOptions;
    acaOptions.minimumBlockSize = 4;
    assemblyOptions.switchToAcaMode(acaOptions);
    AccuracyOptions accuracyOptions;
    accuracyOptions.doubleRegular.setRelativeQuadratureOrder(4);
    accuracyOptions.doubleSingular.setRelativeQuadratureOrder(2);
    NumericalQuadratureStrategy<BFT, RT> quadStrategy(accuracyOptions);

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

    const RT waveNumber = initWaveNumber<RT>();

    BoundaryOperator<BFT, RT> opNonsymmetric =
        modifiedHelmholtz3dSingleLayerBoundaryOperator<BFT, RT, RT>(
            make_shared_from_ref(context),
            make_shared_from_ref(pwiseLinears),
            make_shared_from_ref(pwiseConstants),
            make_shared_from_ref(pwiseLinears),
            waveNumber,
            "", NO_SYMMETRY);
    BoundaryOperator<BFT, RT> opSymmetric =
        modifiedHelmholtz3dSingleLayerBoundaryOperator<BFT, RT, RT>(
            make_shared_from_ref(context),
            make_shared_from_ref(pwiseLinears),
            make_shared_from_ref(pwiseConstants),
            make_shared_from_ref(pwiseLinears),
            waveNumber,
            "", SYMMETRIC);

    arma::Mat<RT> matNonsymmetric = opNonsymmetric.weakForm()->asMatrix();
    arma::Mat<RT> matSymmetric = opSymmetric.weakForm()->asMatrix();

    BOOST_CHECK(check_arrays_are_close<RT>(
                    matNonsymmetric, matSymmetric, 2 * acaOptions.eps));
}

BOOST_AUTO_TEST_CASE_TEMPLATE(aca_of_synthetic_maxwell_single_layer_operator_agrees_with_dense_assembly_in_asymmetric_case,
                              ValueType, complex_result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    GridParameters params;
    params.topology = GridParameters::TRIANGULAR;
    shared_ptr<Grid> grid = GridFactory::importGmshGrid(
        params, "../../meshes/sphere-h-0.4.msh", false /* verbose */);

    RT waveNumber = initWaveNumber<RT>();

    shared_ptr<Space<BFT> > vectorPwiseLinears(
        new RaviartThomas0VectorSpace<BFT>(grid));
    shared_ptr<Space<BFT> > vectorPwiseLinears2(
        new RaviartThomas0VectorSpace<BFT>(grid));

    AccuracyOptions accuracyOptions;
    accuracyOptions.doubleRegular.setRelativeQuadratureOrder(2);
    accuracyOptions.singleRegular.setRelativeQuadratureOrder(2);
    shared_ptr<NumericalQuadratureStrategy<BFT, RT> > quadStrategy(
                new NumericalQuadratureStrategy<BFT, RT>(accuracyOptions));

    AssemblyOptions assemblyOptionsDense;
    assemblyOptionsDense.setVerbosityLevel(VerbosityLevel::LOW);
    shared_ptr<Context<BFT, RT> > contextDense(
        new Context<BFT, RT>(quadStrategy, assemblyOptionsDense));

    BoundaryOperator<BFT, RT> opDense =
            maxwell3dSingleLayerBoundaryOperator<BFT>(
                contextDense,
                vectorPwiseLinears, vectorPwiseLinears, vectorPwiseLinears,
                waveNumber);
    arma::Mat<RT> weakFormDense = opDense.weakForm()->asMatrix();

    AssemblyOptions assemblyOptionsAca;
    assemblyOptionsAca.setVerbosityLevel(VerbosityLevel::LOW);
    AcaOptions acaOptions;
    acaOptions.mode = AcaOptions::LOCAL_ASSEMBLY;
    assemblyOptionsAca.switchToAcaMode(acaOptions);
    shared_ptr<Context<BFT, RT> > contextAca(
        new Context<BFT, RT>(quadStrategy, assemblyOptionsAca));

    // Internal domain different from dualToRange
    BoundaryOperator<BFT, RT> opAca =
            maxwell3dSingleLayerBoundaryOperator<BFT>(
                contextAca,
                vectorPwiseLinears, vectorPwiseLinears, vectorPwiseLinears2,
                waveNumber);
    arma::Mat<RT> weakFormAca = opAca.weakForm()->asMatrix();

    BOOST_CHECK(check_arrays_are_close<ValueType>(
                    weakFormDense, weakFormAca, 2. * acaOptions.eps));
}

BOOST_AUTO_TEST_CASE_TEMPLATE(aca_of_synthetic_modified_helmholtz_hypersingular_operator_agrees_with_dense_assembly_in_symmetric_case,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    GridParameters params;
    params.topology = GridParameters::TRIANGULAR;
    shared_ptr<Grid> grid = GridFactory::importGmshGrid(
        params, "../../meshes/sphere-h-0.4.msh", false /* verbose */);

    RT waveNumber = initWaveNumber<RT>();

    shared_ptr<Space<BFT> > pwiseConstants(
        new PiecewiseConstantScalarSpace<BFT>(grid));
    shared_ptr<Space<BFT> > pwiseLinears(
        new PiecewiseLinearContinuousScalarSpace<BFT>(grid));

    AccuracyOptions accuracyOptions;
    accuracyOptions.doubleRegular.setRelativeQuadratureOrder(2);
    accuracyOptions.singleRegular.setRelativeQuadratureOrder(2);
    shared_ptr<NumericalQuadratureStrategy<BFT, RT> > quadStrategy(
                new NumericalQuadratureStrategy<BFT, RT>(accuracyOptions));

    AssemblyOptions assemblyOptionsDense;
    assemblyOptionsDense.setVerbosityLevel(VerbosityLevel::LOW);
    shared_ptr<Context<BFT, RT> > contextDense(
        new Context<BFT, RT>(quadStrategy, assemblyOptionsDense));

    BoundaryOperator<BFT, RT> opDense =
            modifiedHelmholtz3dHypersingularBoundaryOperator<BFT, RT, RT>(
                contextDense, pwiseLinears, pwiseConstants, pwiseLinears,
                waveNumber);
    arma::Mat<RT> weakFormDense = opDense.weakForm()->asMatrix();

    AssemblyOptions assemblyOptionsAca;
    assemblyOptionsAca.setVerbosityLevel(VerbosityLevel::LOW);
    AcaOptions acaOptions;
    acaOptions.mode = AcaOptions::LOCAL_ASSEMBLY;
    assemblyOptionsAca.switchToAcaMode(acaOptions);
    shared_ptr<Context<BFT, RT> > contextAca(
        new Context<BFT, RT>(quadStrategy, assemblyOptionsAca));

    BoundaryOperator<BFT, RT> opAca =
            modifiedHelmholtz3dHypersingularBoundaryOperator<BFT, RT, RT>(
                contextAca, pwiseLinears, pwiseConstants, pwiseLinears,
                waveNumber);
    arma::Mat<RT> weakFormAca = opAca.weakForm()->asMatrix();

    BOOST_CHECK(check_arrays_are_close<ValueType>(
                    weakFormDense, weakFormAca, 2. * acaOptions.eps));
}

int main()
{
    // Import symbols from namespace Bempp to the global namespace

    using namespace Bempp;

    // Load mesh

    const char* meshFile = "meshes/sphere-h-0.2.msh";
    GridParameters params;
    params.topology = GridParameters::TRIANGULAR;
    shared_ptr<Grid> grid = GridFactory::importGmshGrid(params, meshFile);

    // Initialize the spaces

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

    // Define the quadrature strategy

    AccuracyOptions accuracyOptions;
    // Increase by 2 the order of quadrature rule used to approximate
    // integrals of regular functions on pairs on elements
    accuracyOptions.doubleRegular.setRelativeQuadratureOrder(2);
    // Increase by 2 the order of quadrature rule used to approximate
    // integrals of regular functions on single elements
    accuracyOptions.singleRegular.setRelativeQuadratureOrder(2);
    NumericalQuadratureStrategy<BFT, RT> quadStrategy(accuracyOptions);

    // Specify the assembly method. We want to use ACA

    AssemblyOptions assemblyOptions;
    AcaOptions acaOptions; // Default parameters for ACA
    assemblyOptions.switchToAcaMode(acaOptions);

    // Create the assembly context

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

    // Construct elementary operators

    BoundaryOperator<BFT, RT> slpOp =
            laplace3dSingleLayerBoundaryOperator<BFT, RT>(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseConstants),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseConstants));
    BoundaryOperator<BFT, RT> dlpOp =
            laplace3dDoubleLayerBoundaryOperator<BFT, RT>(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseConstants));
    BoundaryOperator<BFT, RT> idOp =
            identityOperator<BFT, RT>(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseConstants));

    // Form the right-hand side sum

    BoundaryOperator<BFT, RT> rhsOp = -0.5 * idOp + dlpOp;

    // Construct the grid function representing the (input) Dirichlet data

    GridFunction<BFT, RT> dirichletData(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                surfaceNormalIndependentFunction(DirichletData()));

    // Construct the right-hand-side grid function

    GridFunction<BFT, RT> rhs = rhsOp * dirichletData;

    // Initialize the solver

    DefaultIterativeSolver<BFT, RT> solver(slpOp);
    solver.initializeSolver(defaultGmresParameterList(1e-5));

    // Solve the equation

    Solution<BFT, RT> solution = solver.solve(rhs);
    std::cout << solution.solverMessage() << std::endl;

    // Extract the solution in the form of a grid function
    // and export it in VTK format

    const GridFunction<BFT, RT>& solFun = solution.gridFunction();
    solFun.exportToVtk(VtkWriter::CELL_DATA, "Neumann_data", "solution");

    // Compare the numerical and analytical solution on the grid

    // GridFunction<BFT, RT> exactSolFun(
    //             make_shared_from_ref(context),
    //             make_shared_from_ref(pwiseConstants),
    //             make_shared_from_ref(pwiseConstants),
    //             surfaceNormalIndependentFunction(ExactNeumannData()));
    CT absoluteError, relativeError;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    // Physical parameters, general
    const BFT c0 = 0.3;      // speed of light in vacuum [mm/ps]
    BFT refind = 1.4; // refractive index
    BFT alpha = A_Keijzer(refind); // boundary term
    BFT c = c0/refind;       // speed of light in medium [mm/ps]
    BFT freq = 100*1e6; // modulation frequency [Hz]
    BFT omega = 2.0*M_PI * freq*1e-12; // modulation frequency [cycles/ps]

    // Physical parameters, outer region
    BFT mua1 = 0.01; // absorption coefficient
    BFT mus1 = 1.0;  // scattering coefficient
    BFT kappa1 = 1.0/(3.0*(mua1+mus1));   // diffusion coefficient
    RT waveNumber1 = sqrt (RT(mua1/kappa1, omega/(c*kappa1))); // outer region

    // Physical parameters, inner region
    BFT mua2 = 0.02; // absorption coefficient
    BFT mus2 = 0.5;  // scattering coefficient
    BFT kappa2 = 1.0/(3.0*(mua2+mus2));   // diffusion coefficient
    RT waveNumber2 = sqrt (RT(mua2/kappa2, omega/(c*kappa2))); // outer region

    // Create sphere meshes on the fly
    shared_ptr<Grid> grid1 = CreateSphere(25.0, 1.0);
    shared_ptr<Grid> grid2 = CreateSphere(15.0, 1.0);

    // Initialize the spaces

    PiecewiseLinearContinuousScalarSpace<BFT> HplusHalfSpace1(grid1);
    PiecewiseLinearContinuousScalarSpace<BFT> HplusHalfSpace2(grid2);

    // Define some default options.

    AssemblyOptions assemblyOptions;

    // We want to use ACA

    AcaOptions acaOptions; // Default parameters for ACA
    acaOptions.eps = 1e-5;
    assemblyOptions.switchToAcaMode(acaOptions);

    // Define the standard integration factory

    AccuracyOptions accuracyOptions;
    accuracyOptions.doubleRegular.setRelativeQuadratureOrder(2);
    accuracyOptions.singleRegular.setRelativeQuadratureOrder(1);
    NumericalQuadratureStrategy<BFT, RT> quadStrategy(accuracyOptions);
    Context<BFT, RT> context(make_shared_from_ref(quadStrategy), assemblyOptions);

    // We need the single layer, double layer, and the identity operator

    // mesh1 x mesh1
    BoundaryOperator<BFT, RT> slp11 =
            modifiedHelmholtz3dSingleLayerBoundaryOperator<BFT, RT, RT>(
                SHARED(context), SHARED(HplusHalfSpace1),
                SHARED(HplusHalfSpace1), SHARED(HplusHalfSpace1), waveNumber1);
    BoundaryOperator<BFT, RT> dlp11 =
            modifiedHelmholtz3dDoubleLayerBoundaryOperator<BFT, RT, RT>(
                SHARED(context), SHARED(HplusHalfSpace1),
                SHARED(HplusHalfSpace1), SHARED(HplusHalfSpace1), waveNumber1);
    BoundaryOperator<BFT, RT> id11 =
            identityOperator<BFT, RT>(
                SHARED(context), SHARED(HplusHalfSpace1),
                SHARED(HplusHalfSpace1), SHARED(HplusHalfSpace1));

    // mesh2 x mesh2, wavenumber 1
    BoundaryOperator<BFT, RT> slp22_w1 =
            modifiedHelmholtz3dSingleLayerBoundaryOperator<BFT, RT, RT>(
                SHARED(context), SHARED(HplusHalfSpace2),
                SHARED(HplusHalfSpace2), SHARED(HplusHalfSpace2), waveNumber1);
    BoundaryOperator<BFT, RT> dlp22_w1 =
            modifiedHelmholtz3dDoubleLayerBoundaryOperator<BFT, RT, RT>(
                SHARED(context), SHARED(HplusHalfSpace2),
                SHARED(HplusHalfSpace2), SHARED(HplusHalfSpace2), waveNumber1);
    BoundaryOperator<BFT, RT> id22 =
            identityOperator<BFT, RT>(
                SHARED(context), SHARED(HplusHalfSpace2),
                SHARED(HplusHalfSpace2), SHARED(HplusHalfSpace2));

    // mesh2 x mesh2, wavenumber 2
    BoundaryOperator<BFT, RT> slp22_w2 =
            modifiedHelmholtz3dSingleLayerBoundaryOperator<BFT, RT, RT>(
                SHARED(context), SHARED(HplusHalfSpace2),
                SHARED(HplusHalfSpace2), SHARED(HplusHalfSpace2), waveNumber2);
    BoundaryOperator<BFT, RT> dlp22_w2 =
            modifiedHelmholtz3dDoubleLayerBoundaryOperator<BFT, RT, RT>(
                SHARED(context), SHARED(HplusHalfSpace2),
                SHARED(HplusHalfSpace2), SHARED(HplusHalfSpace2), waveNumber2);

    // mesh1 x mesh2
    BoundaryOperator<BFT, RT> slp12 =
            modifiedHelmholtz3dSingleLayerBoundaryOperator<BFT, RT, RT>(
                SHARED(context), SHARED(HplusHalfSpace2),
                SHARED(HplusHalfSpace1), SHARED(HplusHalfSpace1), waveNumber1);
    BoundaryOperator<BFT, RT> dlp12 =
            modifiedHelmholtz3dDoubleLayerBoundaryOperator<BFT, RT, RT>(
                SHARED(context), SHARED(HplusHalfSpace2),
                SHARED(HplusHalfSpace1), SHARED(HplusHalfSpace1), waveNumber1);

    // mesh2 x mesh1
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    // Load mesh

    if (argc != 2) {
        std::cout << "Solve a Dirichlet problem for the Laplace equation.\n"
                     "Usage: " << argv[0] << " <mesh_file>" << std::endl;
        return 1;
    }
    shared_ptr<Grid> grid = loadTriangularMeshFromFile(argv[1]);

    //std::cout << grid->gridTopology() << std::endl;

    // Initialize the spaces

    PiecewiseLinearContinuousScalarSpace<BFT> HplusHalfSpace(grid);
    PiecewiseConstantScalarSpace<BFT> HminusHalfSpace(grid);

    // Define some default options.

    AssemblyOptions assemblyOptions;

    // We want to use ACA

    AcaOptions acaOptions; // Default parameters for ACA
    acaOptions.eps = 1e-5;
    assemblyOptions.switchToAcaMode(acaOptions);

    // Define the standard integration factory

    AccuracyOptions accuracyOptions;
    accuracyOptions.doubleRegular.setRelativeQuadratureOrder(1);
    NumericalQuadratureStrategy<BFT, RT> quadStrategy(accuracyOptions);

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

    // We need the single layer, double layer, and the identity operator
    BoundaryOperator<BFT, RT> slpOp = laplace3dSingleLayerBoundaryOperator<BFT, RT >(
                make_shared_from_ref(context),
                make_shared_from_ref(HminusHalfSpace),
                make_shared_from_ref(HplusHalfSpace),
                make_shared_from_ref(HminusHalfSpace),
                "SLP");
    BoundaryOperator<BFT, RT> dlpOp = laplace3dDoubleLayerBoundaryOperator<BFT, RT >(
                make_shared_from_ref(context),
                make_shared_from_ref(HplusHalfSpace),
                make_shared_from_ref(HplusHalfSpace),
                make_shared_from_ref(HminusHalfSpace),
                "DLP");
    BoundaryOperator<BFT, RT> id = identityOperator<BFT, RT>(
                make_shared_from_ref(context),
                make_shared_from_ref(HplusHalfSpace),
                make_shared_from_ref(HplusHalfSpace),
                make_shared_from_ref(HminusHalfSpace),
                "I");

    // Form the right-hand side sum

    BoundaryOperator<BFT, RT> rhsOp = -0.5 * id + dlpOp;

    // We also want a grid function

    GridFunction<BFT, RT> u(
                make_shared_from_ref(context),
                make_shared_from_ref(HplusHalfSpace),
                make_shared_from_ref(HplusHalfSpace), // is this the right choice?
                surfaceNormalIndependentFunction(MyFunctor()));

    // Assemble the rhs

    std::cout << "Assemble rhs" << std::endl;

    GridFunction<BFT, RT> rhs = rhsOp * u;

    // Initialize the solver

    std::cout << "Initialize solver" << std::endl;

#ifdef WITH_TRILINOS
    DefaultIterativeSolver<BFT, RT> solver(slpOp,ConvergenceTestMode::TEST_CONVERGENCE_IN_DUAL_TO_RANGE);
    solver.initializeSolver(defaultGmresParameterList(1e-5));
    Solution<BFT, RT> solution = solver.solve(rhs);
    std::cout << solution.solverMessage() << std::endl;
#else
    DefaultDirectSolver<BFT, RT> solver(slp, rhs);
    solver.solve();
#endif

    // Extract the solution

    const GridFunction<BFT, RT>& solFun = solution.gridFunction();

    // Write out as VTK

    solFun.exportToVtk(VtkWriter::CELL_DATA, "Neumann_data",
                       "calculated_neumann_data_cell");
    solFun.exportToVtk(VtkWriter::VERTEX_DATA, "Neumann_data",
                       "calculated_neumann_data_vertex");

    // Uncomment the block below if you are solving the problem on a sphere and
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    // Process command-line args

    if (argc < 7 || argc % 2 != 1) {
        std::cout << "Solve a Maxwell Dirichlet problem in an exterior domain.\n"
                  << "Usage: " << argv[0]
                  << " <mesh_file> <n_threads> <aca_eps> <solver_tol>"
                  << " <singular_order_increment>"
                  << " [<regular_order_increment_1> <min_relative_distance_1>]"
                  << " [<regular_order_increment_2> <min_relative_distance_2>]"
                  << " [...] <regular_order_increment_n>"
                  << std::endl;
        return 1;
    }
    int maxThreadCount = atoi(argv[2]);
    double acaEps = atof(argv[3]);
    double solverTol = atof(argv[4]);
    int singOrderIncrement = atoi(argv[5]);
    if (acaEps > 1. || acaEps < 0.) {
        std::cout << "Invalid aca_eps: " << acaEps << std::endl;
        return 1;
    }
    if (solverTol > 1. || solverTol < 0.) {
        std::cout << "Invalid solver_tol: " << solverTol << std::endl;
        return 1;
    }

    AccuracyOptionsEx accuracyOptions;
    std::vector<double> maxNormDists;
    std::vector<int> orderIncrements;
    for (int i = 6; i < argc - 1; i += 2) {
        orderIncrements.push_back(atoi(argv[i]));
        maxNormDists.push_back(atof(argv[i + 1]));
    }
    orderIncrements.push_back(atoi(argv[argc - 1]));
    accuracyOptions.setDoubleRegular(maxNormDists, orderIncrements);
    accuracyOptions.setDoubleSingular(singOrderIncrement);
    accuracyOptions.setSingleRegular(2);

    // Load mesh

    GridParameters params;
    params.topology = GridParameters::TRIANGULAR;
    shared_ptr<Grid> grid = GridFactory::importGmshGrid(params, argv[1]);

    // Initialize the space

    RaviartThomas0VectorSpace<BFT> HdivSpace(grid);

    // Set assembly mode and options

    AssemblyOptions assemblyOptions;
    assemblyOptions.setMaxThreadCount(maxThreadCount);
    if (acaEps > 0.) {
        AcaOptions acaOptions;
        acaOptions.eps = acaEps;
        assemblyOptions.switchToAcaMode(acaOptions);
    }

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

    // Construct operators

    BoundaryOperator<BFT, RT> slpOp = maxwell3dSingleLayerBoundaryOperator<BFT>(
                make_shared_from_ref(context),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                k,
                "SLP");
    BoundaryOperator<BFT, RT> dlpOp = maxwell3dDoubleLayerBoundaryOperator<BFT>(
                make_shared_from_ref(context),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                k,
                "DLP");
    BoundaryOperator<BFT, RT> idOp = maxwell3dIdentityOperator<BFT, RT>(
                make_shared_from_ref(context),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                "Id");

    // Construct a grid function representing the Dirichlet data

    GridFunction<BFT, RT> dirichletData(
                make_shared_from_ref(context),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                surfaceNormalDependentFunction(DirichletData()));

    dirichletData.exportToVtk(VtkWriter::CELL_DATA, "Dirichlet_data",
                              "input_dirichlet_data_cell");
    dirichletData.exportToVtk(VtkWriter::VERTEX_DATA, "Dirichlet_data",
                              "input_dirichlet_data_vertex");

    // Construct a grid function representing the right-hand side
//.........这里部分代码省略.........