Skip to content

Commit

Permalink
Sorting out boundary conditions (#28)
Browse files Browse the repository at this point in the history 
* 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.
  • Loading branch information
@dreamer2368
dreamer2368 authored Mar 4, 2024
1 parent ed5b9a8 commit 6fa043d
Show file tree
Hide file tree
Showing 10 changed files with 131 additions and 80 deletions.
3 changes: 0 additions & 3 deletions include/linelast_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,9 +52,6 @@ class LinElastSolver : public MultiBlockSolver
// Initial positions
VectorFunctionCoefficient *init_x = NULL;

// Boundary condition types
Array<LinElastProblem::BoundaryType> type_idx;

public:
LinElastSolver();

Expand Down
4 changes: 3 additions & 1 deletion include/multiblock_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,7 @@ friend class ParameterizedProblem;
int max_bdr_attr;
int numBdr;
Array<Array<int> *> bdr_markers;
Array<BoundaryType> bdr_type; // Boundary condition types of (numBdr) array

// MFEM solver options
bool use_amg;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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<GridFunction *> fom_sols, Array<GridFunction *> rom_sols, Vector &error);
Expand Down
20 changes: 9 additions & 11 deletions include/parameterized_problem.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -128,7 +136,7 @@ friend class RandomSampleGenerator;
function_factory::GeneralVectorFunction *vector_rhs_ptr = NULL;
Array<function_factory::GeneralVectorFunction *> vector_bdr_ptr;
Array<int> battr;
Array<int> bdr_type; // abstract boundary type
Array<BoundaryType> bdr_type; // abstract boundary type

Array<function_factory::GeneralScalarFunction *> general_scalar_ptr;
Array<function_factory::GeneralVectorFunction *> general_vector_ptr;
Expand All @@ -145,10 +153,6 @@ class PoissonProblem : public ParameterizedProblem
{
friend class PoissonSolver;

protected:
enum BoundaryType
{ ZERO, DIRICHLET, NEUMANN, NUM_BDR_TYPE };

public:
virtual ~PoissonProblem() {};
};
Expand Down Expand Up @@ -183,10 +187,6 @@ class StokesProblem : public ParameterizedProblem
{
friend class StokesSolver;

protected:
enum BoundaryType
{ ZERO, DIRICHLET, NEUMANN, NUM_BDR_TYPE };

public:
virtual ~StokesProblem() {};
};
Expand Down Expand Up @@ -222,8 +222,6 @@ class LinElastProblem : public ParameterizedProblem
friend class LinElastSolver;

protected:
enum BoundaryType
{ ZERO, DIRICHLET, NEUMANN, NUM_BDR_TYPE };
Vector lambda, mu;

public:
Expand Down
59 changes: 31 additions & 28 deletions src/linelast_solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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()
Expand Down Expand Up @@ -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;
}
}
}
Expand Down Expand Up @@ -427,14 +428,20 @@ 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]));
}
}
}
}

void LinElastSolver::SetParameterizedProblem(ParameterizedProblem *problem)
{
/* set up boundary types */
MultiBlockSolver::SetParameterizedProblem(problem);

// Set materials
lambda_c.SetSize(numSub);
lambda_c = NULL;
Expand All @@ -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
Expand Down
31 changes: 31 additions & 0 deletions src/multiblock_solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<FiniteElementSpace *> &comp_fes)
Expand Down Expand Up @@ -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++)
Expand Down Expand Up @@ -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;
Expand Down
Loading

0 comments on commit 6fa043d

Please sign in to comment.