From 6fa043d1cb13316177155e1f4a55931334b1dcc8 Mon Sep 17 00:00:00 2001 From: "Kevin\" Seung Whan Chung" Date: Mon, 4 Mar 2024 11:29:45 -0800 Subject: [PATCH] Sorting out boundary conditions (#28) * BoundaryType belongs to ParameterizedProblem * bdr type available for all MultiBlockSolver class. * MultiBlockSolver::SetBdrType * change on bc handling * BoundaryType::ZERO is reserved for zero dirichlet bc. * BoundaryType moved up as independent enumeration * allow zero neumann condition for linear elasticity. --- include/linelast_solver.hpp | 3 -- include/multiblock_solver.hpp | 4 ++- include/parameterized_problem.hpp | 20 +++++------ src/linelast_solver.cpp | 59 ++++++++++++++++--------------- src/multiblock_solver.cpp | 31 ++++++++++++++++ src/parameterized_problem.cpp | 57 +++++++++++++++-------------- src/poisson_solver.cpp | 10 ++++-- src/steady_ns_solver.cpp | 6 ++-- src/stokes_solver.cpp | 17 ++++++--- test/mms_suite.cpp | 4 +++ 10 files changed, 131 insertions(+), 80 deletions(-) diff --git a/include/linelast_solver.hpp b/include/linelast_solver.hpp index 819dc60e..0ae4b356 100644 --- a/include/linelast_solver.hpp +++ b/include/linelast_solver.hpp @@ -52,9 +52,6 @@ class LinElastSolver : public MultiBlockSolver // Initial positions VectorFunctionCoefficient *init_x = NULL; - // Boundary condition types - Array type_idx; - public: LinElastSolver(); diff --git a/include/multiblock_solver.hpp b/include/multiblock_solver.hpp index a9006c86..1641f2f4 100644 --- a/include/multiblock_solver.hpp +++ b/include/multiblock_solver.hpp @@ -68,6 +68,7 @@ friend class ParameterizedProblem; int max_bdr_attr; int numBdr; Array *> bdr_markers; + Array bdr_type; // Boundary condition types of (numBdr) array // MFEM solver options bool use_amg; @@ -153,6 +154,7 @@ friend class ParameterizedProblem; virtual void BuildRHSOperators() = 0; virtual void BuildDomainOperators() = 0; + void SetBdrType(const BoundaryType type, const int &global_battr_idx=-1); virtual bool BCExistsOnBdr(const int &global_battr_idx) = 0; virtual void SetupBCOperators() = 0; virtual void SetupRHSBCOperators() = 0; @@ -235,7 +237,7 @@ friend class ParameterizedProblem; virtual void SaveBasisVisualization() { rom_handler->SaveBasisVisualization(fes, var_names); } - virtual void SetParameterizedProblem(ParameterizedProblem *problem) = 0; + virtual void SetParameterizedProblem(ParameterizedProblem *problem); void ComputeSubdomainErrorAndNorm(GridFunction *fom_sol, GridFunction *rom_sol, double &error, double &norm); void ComputeRelativeError(Array fom_sols, Array rom_sols, Vector &error); diff --git a/include/parameterized_problem.hpp b/include/parameterized_problem.hpp index bc210eb4..e815a3e2 100644 --- a/include/parameterized_problem.hpp +++ b/include/parameterized_problem.hpp @@ -94,6 +94,14 @@ void tip_force(const Vector &x, Vector &u); } +enum BoundaryType +{ + ZERO, + DIRICHLET, + NEUMANN, + NUM_BDR_TYPE +}; + class ParameterizedProblem { friend class SampleGenerator; @@ -128,7 +136,7 @@ friend class RandomSampleGenerator; function_factory::GeneralVectorFunction *vector_rhs_ptr = NULL; Array vector_bdr_ptr; Array battr; - Array bdr_type; // abstract boundary type + Array bdr_type; // abstract boundary type Array general_scalar_ptr; Array general_vector_ptr; @@ -145,10 +153,6 @@ class PoissonProblem : public ParameterizedProblem { friend class PoissonSolver; -protected: - enum BoundaryType - { ZERO, DIRICHLET, NEUMANN, NUM_BDR_TYPE }; - public: virtual ~PoissonProblem() {}; }; @@ -183,10 +187,6 @@ class StokesProblem : public ParameterizedProblem { friend class StokesSolver; -protected: - enum BoundaryType - { ZERO, DIRICHLET, NEUMANN, NUM_BDR_TYPE }; - public: virtual ~StokesProblem() {}; }; @@ -222,8 +222,6 @@ class LinElastProblem : public ParameterizedProblem friend class LinElastSolver; protected: - enum BoundaryType - { ZERO, DIRICHLET, NEUMANN, NUM_BDR_TYPE }; Vector lambda, mu; public: diff --git a/src/linelast_solver.cpp b/src/linelast_solver.cpp index 821695a3..d7e7259c 100644 --- a/src/linelast_solver.cpp +++ b/src/linelast_solver.cpp @@ -75,9 +75,6 @@ void LinElastSolver::SetupBCVariables() bdr_coeffs.SetSize(numBdr); bdr_coeffs = NULL; - type_idx.SetSize(numBdr); - type_idx = LinElastProblem::BoundaryType::DIRICHLET; - lambda_c.SetSize(numSub); lambda_c = NULL; @@ -144,14 +141,12 @@ void LinElastSolver::BuildOperators() BuildRHSOperators(); BuildDomainOperators(); } + bool LinElastSolver::BCExistsOnBdr(const int &global_battr_idx) { assert((global_battr_idx >= 0) && (global_battr_idx < global_bdr_attributes.Size())); assert(bdr_coeffs.Size() == global_bdr_attributes.Size()); - if(type_idx[global_battr_idx] == LinElastProblem::BoundaryType::NEUMANN) - return false; - else - return (bdr_coeffs[global_battr_idx]); + return (bdr_coeffs[global_battr_idx]); } void LinElastSolver::BuildRHSOperators() @@ -179,17 +174,23 @@ void LinElastSolver::SetupRHSBCOperators() if (idx < 0) continue; - if (type_idx[b] == LinElastProblem::BoundaryType::NEUMANN) - { - bs[m]->AddBdrFaceIntegrator(new VectorBoundaryLFIntegrator(*bdr_coeffs[b]), *bdr_markers[b]); - } - if (!BCExistsOnBdr(b)) continue; - if (type_idx[b] == LinElastProblem::BoundaryType::DIRICHLET) // Default behavior - { - bs[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator(*bdr_coeffs[b], *lambda_c[m], *mu_c[m], alpha, kappa), *bdr_markers[b]); + switch (bdr_type[b]) + { + case BoundaryType::DIRICHLET: + bs[m]->AddBdrFaceIntegrator(new DGElasticityDirichletLFIntegrator( + *bdr_coeffs[b], *lambda_c[m], *mu_c[m], alpha, kappa), *bdr_markers[b]); + break; + case BoundaryType::NEUMANN: + bs[m]->AddBdrFaceIntegrator(new VectorBoundaryLFIntegrator(*bdr_coeffs[b]), *bdr_markers[b]); + break; + default: + printf("LinElastSolver::SetupRHSBCOperators - "); + printf("boundary attribute %d has a non-zero function, but does not have boundary type. will not be enforced.", + global_bdr_attributes[b]); + break; } } } @@ -427,7 +428,10 @@ void LinElastSolver::SetupDomainBCOperators() continue; if (!BCExistsOnBdr(b)) continue; - as[m]->AddBdrFaceIntegrator(new DGElasticityIntegrator(*(lambda_c[m]), *(mu_c[m]), alpha, kappa), *(bdr_markers[b])); + + if (bdr_type[b] == BoundaryType::DIRICHLET) + as[m]->AddBdrFaceIntegrator(new DGElasticityIntegrator( + *(lambda_c[m]), *(mu_c[m]), alpha, kappa), *(bdr_markers[b])); } } } @@ -435,6 +439,9 @@ void LinElastSolver::SetupDomainBCOperators() void LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem) { + /* set up boundary types */ + MultiBlockSolver::SetParameterizedProblem(problem); + // Set materials lambda_c.SetSize(numSub); lambda_c = NULL; @@ -454,19 +461,15 @@ void LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem) } // Set BCs, the switch on BC type is done inside SetupRHSBCOperators - for (int b = 0; b < global_bdr_attributes.Size(); b++) + for (int b = 0; b < problem->battr.Size(); b++) { - int ti = problem->battr.Find(global_bdr_attributes[b]); - - if (ti >=0) - { - type_idx[b] = (LinElastProblem::BoundaryType)problem->bdr_type[ti]; - if (problem->bdr_type[ti] == LinElastProblem::BoundaryType::DIRICHLET || problem->bdr_type[ti] == LinElastProblem::BoundaryType::NEUMANN) - { - assert(problem->vector_bdr_ptr[ti]); - AddBCFunction(*(problem->vector_bdr_ptr[ti]), problem->battr[ti]); - } - } + /* Dirichlet bc requires a function specified, even for zero. */ + if (problem->bdr_type[b] == BoundaryType::DIRICHLET) + assert(problem->vector_bdr_ptr[b]); + + /* Neumann bc does not require a function specified for zero */ + if (problem->vector_bdr_ptr[b]) + AddBCFunction(*(problem->vector_bdr_ptr[b]), problem->battr[b]); } // Set RHS diff --git a/src/multiblock_solver.cpp b/src/multiblock_solver.cpp index f3bdf1d4..0d4ae09d 100644 --- a/src/multiblock_solver.cpp +++ b/src/multiblock_solver.cpp @@ -218,6 +218,21 @@ void MultiBlockSolver::SetupBCVariables() (*bdr_markers[k]) = 0; (*bdr_markers[k])[bdr_attr-1] = 1; } + + bdr_type.SetSize(global_bdr_attributes.Size()); + bdr_type = BoundaryType::NUM_BDR_TYPE; +} + +void MultiBlockSolver::SetBdrType(const BoundaryType type, const int &global_battr_idx) +{ + assert(bdr_type.Size() == global_bdr_attributes.Size()); + if (global_battr_idx < 0) + bdr_type = type; + else + { + assert(global_battr_idx < global_bdr_attributes.Size()); + bdr_type[global_battr_idx] = type; + } } void MultiBlockSolver::GetComponentFESpaces(Array &comp_fes) @@ -542,6 +557,10 @@ void MultiBlockSolver::AssembleROM() if (global_idx < 0) continue; if (!BCExistsOnBdr(global_idx)) continue; + /* we assume only Neumann condition would not add an operator. */ + if (bdr_type[global_idx] == BoundaryType::NEUMANN) + continue; + MatrixBlocks *bdr_mat = (*bdr_mats[c_type])[b]; AddToBlockMatrix(midx, midx, *bdr_mat, *romMat); } // for (int b = 0; b < bdr_c2g->Size(); b++) @@ -685,6 +704,18 @@ void MultiBlockSolver::SaveVisualization() } }; +void MultiBlockSolver::SetParameterizedProblem(ParameterizedProblem *problem) +{ + assert(bdr_type.Size() == global_bdr_attributes.Size()); + for (int b = 0; b < global_bdr_attributes.Size(); b++) + { + int idx = problem->battr.Find(global_bdr_attributes[b]); + if (idx < 0) continue; + + bdr_type[b] = problem->bdr_type[idx]; + } +} + void MultiBlockSolver::SaveSolution(std::string filename) { if (!save_sol) return; diff --git a/src/parameterized_problem.cpp b/src/parameterized_problem.cpp index e24e5a67..a40119e2 100644 --- a/src/parameterized_problem.cpp +++ b/src/parameterized_problem.cpp @@ -206,7 +206,7 @@ ParameterizedProblem::ParameterizedProblem() : problem_name(config.GetRequiredOption("parameterized_problem/name")) { battr.SetSize(1); battr = -1; - bdr_type.SetSize(1); bdr_type = -1; + bdr_type.SetSize(1); bdr_type = BoundaryType::NUM_BDR_TYPE; scalar_bdr_ptr.SetSize(1); vector_bdr_ptr.SetSize(1); @@ -320,7 +320,7 @@ Poisson0::Poisson0() { param_num = 2; battr = -1; - bdr_type = PoissonProblem::ZERO; + bdr_type = BoundaryType::ZERO; scalar_bdr_ptr.SetSize(1); vector_bdr_ptr.SetSize(1); @@ -351,7 +351,7 @@ PoissonComponent::PoissonComponent() // k (max 3) + offset (1) + bdr_k (max 3) + bdr_offset(1) + bdr_idx(1) param_num = 9; battr = -1; - bdr_type = PoissonProblem::DIRICHLET; + bdr_type = BoundaryType::DIRICHLET; // pointer to static function. scalar_rhs_ptr = &(function_factory::poisson_component::rhs); @@ -417,7 +417,7 @@ PoissonSpiral::PoissonSpiral() { param_num = 4; battr = -1; - bdr_type = PoissonProblem::ZERO; + bdr_type = BoundaryType::ZERO; // pointer to static function. scalar_bdr_ptr = NULL; @@ -454,8 +454,8 @@ StokesChannel::StokesChannel() battr[2] = 4; battr[3] = 5; bdr_type.SetSize(4); - bdr_type = StokesProblem::ZERO; - bdr_type[2] = StokesProblem::DIRICHLET; + bdr_type = BoundaryType::ZERO; + bdr_type[2] = BoundaryType::DIRICHLET; // pointer to static function. vector_bdr_ptr.SetSize(4); @@ -489,8 +489,8 @@ StokesComponent::StokesComponent() for (int b = 0; b < 5; b++) battr[b] = b+1; bdr_type.SetSize(5); - bdr_type = StokesProblem::DIRICHLET; - bdr_type[4] = StokesProblem::ZERO; + bdr_type = BoundaryType::DIRICHLET; + bdr_type[4] = BoundaryType::ZERO; // pointer to static function. vector_bdr_ptr.SetSize(5); @@ -564,24 +564,24 @@ void StokesFlowPastArray::SetBattr() { if ((*u0)[0] > 0.0) { - bdr_type[3] = StokesProblem::DIRICHLET; - bdr_type[1] = StokesProblem::NEUMANN; + bdr_type[3] = BoundaryType::DIRICHLET; + bdr_type[1] = BoundaryType::NEUMANN; } else { - bdr_type[1] = StokesProblem::DIRICHLET; - bdr_type[3] = StokesProblem::NEUMANN; + bdr_type[1] = BoundaryType::DIRICHLET; + bdr_type[3] = BoundaryType::NEUMANN; } if ((*u0)[1] > 0.0) { - bdr_type[0] = StokesProblem::DIRICHLET; - bdr_type[2] = StokesProblem::NEUMANN; + bdr_type[0] = BoundaryType::DIRICHLET; + bdr_type[2] = BoundaryType::NEUMANN; } else { - bdr_type[2] = StokesProblem::DIRICHLET; - bdr_type[0] = StokesProblem::NEUMANN; + bdr_type[2] = BoundaryType::DIRICHLET; + bdr_type[0] = BoundaryType::NEUMANN; } } @@ -598,13 +598,13 @@ LinElastDisp::LinElastDisp() vector_bdr_ptr.SetSize(3); for (size_t i = 0; i < vector_bdr_ptr.Size(); i++) { - bdr_type[i] = LinElastProblem::DIRICHLET; - battr[i] = i+1; - vector_bdr_ptr[i] = &(function_factory::linelast_disp::init_disp); + bdr_type[i] = BoundaryType::DIRICHLET; + battr[i] = i+1; + vector_bdr_ptr[i] = &(function_factory::linelast_disp::init_disp); } battr[2] = 3; - bdr_type[2] = LinElastProblem::ZERO; + bdr_type[2] = BoundaryType::ZERO; vector_bdr_ptr[2] = NULL; // Set materials @@ -641,12 +641,13 @@ LinElastDispLCantilever::LinElastDispLCantilever() for (size_t i = 0; i < 2; i++) { battr[i] = i+1; - bdr_type[i] = LinElastProblem::DIRICHLET; + bdr_type[i] = BoundaryType::DIRICHLET; vector_bdr_ptr[i] = &(function_factory::linelast_disp::init_disp_lcantilever); } + /* homogeneous Neumann bc */ battr[2] = 3; - bdr_type[2] = LinElastProblem::ZERO; + bdr_type[2] = BoundaryType::NEUMANN; vector_bdr_ptr[2] = NULL; // Set materials @@ -682,12 +683,13 @@ LinElastDispLattice::LinElastDispLattice() for (size_t i = 0; i < 2; i++) { battr[i] = i+1; - bdr_type[i] = LinElastProblem::DIRICHLET; + bdr_type[i] = BoundaryType::DIRICHLET; vector_bdr_ptr[i] = &(function_factory::linelast_disp::init_disp); } + /* homogeneous Neumann bc */ battr[2] = 3; - bdr_type[2] = LinElastProblem::ZERO; + bdr_type[2] = BoundaryType::NEUMANN; vector_bdr_ptr[2] = NULL; // Set materials @@ -725,16 +727,17 @@ LinElastForceCantilever::LinElastForceCantilever() // Fixed end of cantilever battr[0] = 1; - bdr_type[0] = LinElastProblem::DIRICHLET; + bdr_type[0] = BoundaryType::DIRICHLET; vector_bdr_ptr[0] = &(function_factory::linelast_disp::init_disp); // Free end of cantilever battr[1] = 2; - bdr_type[1] = LinElastProblem::NEUMANN; + bdr_type[1] = BoundaryType::NEUMANN; vector_bdr_ptr[1] = &(function_factory::linelast_force::tip_force); + /* homogeneous Neumann bc */ battr[2] = 3; - bdr_type[2] = LinElastProblem::ZERO; + bdr_type[2] = BoundaryType::NEUMANN; vector_bdr_ptr[2] = NULL; // Set materials diff --git a/src/poisson_solver.cpp b/src/poisson_solver.cpp index 7922be6d..297a9a42 100644 --- a/src/poisson_solver.cpp +++ b/src/poisson_solver.cpp @@ -225,6 +225,8 @@ void PoissonSolver::SetupRHSBCOperators() int idx = meshes[m]->bdr_attributes.Find(global_bdr_attributes[b]); if (idx < 0) continue; if (!BCExistsOnBdr(b)) continue; + if (bdr_type[b] == BoundaryType::NEUMANN) + continue; bs[m]->AddBdrFaceIntegrator(new DGDirichletLFIntegrator(*bdr_coeffs[b], sigma, kappa), *bdr_markers[b]); } @@ -245,6 +247,8 @@ void PoissonSolver::SetupDomainBCOperators() int idx = meshes[m]->bdr_attributes.Find(global_bdr_attributes[b]); if (idx < 0) continue; if (!BCExistsOnBdr(b)) continue; + if (bdr_type[b] == BoundaryType::NEUMANN) + continue; as[m]->AddBdrFaceIntegrator(new DGDiffusionIntegrator(sigma, kappa), *bdr_markers[b]); } @@ -565,15 +569,15 @@ void PoissonSolver::SetParameterizedProblem(ParameterizedProblem *problem) { switch (problem->bdr_type[b]) { - case PoissonProblem::BoundaryType::DIRICHLET: + case BoundaryType::DIRICHLET: { assert(problem->scalar_bdr_ptr[b]); AddBCFunction(*(problem->scalar_bdr_ptr[b]), problem->battr[b]); break; } - case PoissonProblem::BoundaryType::NEUMANN: break; + case BoundaryType::NEUMANN: break; default: - case PoissonProblem::BoundaryType::ZERO: + case BoundaryType::ZERO: { AddBCFunction(0.0, problem->battr[b]); break; } } } diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index 56f09242..7282fab0 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -394,8 +394,10 @@ void SteadyNSSolver::SetupRHSBCOperators() { int idx = meshes[m]->bdr_attributes.Find(global_bdr_attributes[b]); if (idx < 0) continue; - // TODO: Non-homogeneous Neumann stress bc if (!BCExistsOnBdr(b)) continue; + // TODO: Non-homogeneous Neumann stress bc + if (bdr_type[b] == BoundaryType::NEUMANN) + continue; fs[m]->AddBdrFaceIntegrator(new DGBdrTemamLFIntegrator(*ud_coeffs[b], minus_half_zeta), *bdr_markers[b]); } @@ -418,7 +420,7 @@ void SteadyNSSolver::SetupDomainBCOperators() if (idx < 0) continue; // homogeneous Neumann boundary condition - if (!BCExistsOnBdr(b)) + if (!BCExistsOnBdr(b) && (bdr_type[b] == BoundaryType::NEUMANN)) hs[m]->AddBdrFaceIntegrator(new DGTemamFluxIntegrator(*zeta_coeff), *bdr_markers[b]); } } diff --git a/src/stokes_solver.cpp b/src/stokes_solver.cpp index e9772c0b..65b4cb50 100644 --- a/src/stokes_solver.cpp +++ b/src/stokes_solver.cpp @@ -324,9 +324,12 @@ void StokesSolver::SetupRHSBCOperators() { int idx = meshes[m]->bdr_attributes.Find(global_bdr_attributes[b]); if (idx < 0) continue; - // TODO: Non-homogeneous Neumann stress bc if (!BCExistsOnBdr(b)) continue; + // TODO: Non-homogeneous Neumann stress bc + if (bdr_type[b] == BoundaryType::NEUMANN) + continue; + fs[m]->AddBdrFaceIntegrator(new DGVectorDirichletLFIntegrator(*ud_coeffs[b], *nu_coeff, sigma, kappa), *bdr_markers[b]); // TODO: Non-homogeneous Neumann stress bc @@ -356,6 +359,10 @@ void StokesSolver::SetupDomainBCOperators() if (idx < 0) continue; if (!BCExistsOnBdr(b)) continue; + // TODO: Non-homogeneous Neumann stress bc + if (bdr_type[b] == BoundaryType::NEUMANN) + continue; + ms[m]->AddBdrFaceIntegrator(new DGVectorDiffusionIntegrator(*nu_coeff, sigma, kappa), *bdr_markers[b]); bs[m]->AddBdrFaceIntegrator(new DGNormalFluxIntegrator, *bdr_markers[b]); } @@ -1093,15 +1100,15 @@ void StokesSolver::SetParameterizedProblem(ParameterizedProblem *problem) { switch (problem->bdr_type[b]) { - case StokesProblem::BoundaryType::DIRICHLET: + case BoundaryType::DIRICHLET: { assert(problem->vector_bdr_ptr[b]); AddBCFunction(*(problem->vector_bdr_ptr[b]), problem->battr[b]); break; } - case StokesProblem::BoundaryType::NEUMANN: break; + case BoundaryType::NEUMANN: break; default: - case StokesProblem::BoundaryType::ZERO: + case BoundaryType::ZERO: { AddBCFunction(zero, problem->battr[b]); break; } } } @@ -1120,7 +1127,7 @@ void StokesSolver::SetParameterizedProblem(ParameterizedProblem *problem) nz_dbcs = true; for (int b = 0; b < problem->battr.Size(); b++) { - if (problem->bdr_type[b] == StokesProblem::BoundaryType::ZERO) + if (problem->bdr_type[b] == BoundaryType::ZERO) { if (problem->battr[b] == -1) { nz_dbcs = false; break; } diff --git a/test/mms_suite.cpp b/test/mms_suite.cpp index c0dd6e87..2d2ee418 100644 --- a/test/mms_suite.cpp +++ b/test/mms_suite.cpp @@ -40,6 +40,7 @@ PoissonSolver *SolveWithRefinement(const int num_refinement) test->InitVisualization(); test->AddBCFunction(ExactSolution); + test->SetBdrType(BoundaryType::DIRICHLET); test->AddRHSFunction(ExactRHS); test->BuildOperators(); @@ -178,6 +179,7 @@ StokesSolver *SolveWithRefinement(const int num_refinement) test->InitVisualization(); test->AddBCFunction(uFun_ex); + test->SetBdrType(BoundaryType::DIRICHLET); test->AddRHSFunction(fFun); // NOTE(kevin): uFun_ex already satisfies zero divergence. // no need to set complementary flux. @@ -357,6 +359,7 @@ SteadyNSSolver *SolveWithRefinement(const int num_refinement) test->InitVisualization(); test->AddBCFunction(uFun_ex); + test->SetBdrType(BoundaryType::DIRICHLET); test->AddRHSFunction(fFun); // NOTE(kevin): uFun_ex already satisfies zero divergence. // no need to set complementary flux. @@ -516,6 +519,7 @@ namespace linelast test->AddBCFunction(ExactSolution, 1); test->AddBCFunction(ExactSolution, 2); test->AddBCFunction(ExactSolution, 3); + test->SetBdrType(BoundaryType::DIRICHLET); test->AddRHSFunction(ExactRHS); test->BuildOperators();