Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Sorting out boundary conditions #28

Merged
merged 7 commits into from
Mar 4, 2024
Merged

Sorting out boundary conditions #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
0a6d19a
BoundaryType belongs to ParameterizedProblem
dreamer2368 Mar 2, 2024
5da4723
bdr type available for all MultiBlockSolver class.
dreamer2368 Mar 2, 2024
fb6ded5
MultiBlockSolver::SetBdrType
dreamer2368 Mar 2, 2024
18b45ff
change on bc handling
dreamer2368 Mar 2, 2024
1d70e11
BoundaryType::ZERO is reserved for zero dirichlet bc.
dreamer2368 Mar 2, 2024
ee877a4
BoundaryType moved up as independent enumeration
dreamer2368 Mar 2, 2024
2f7ee72
allow zero neumann condition for linear elasticity.
dreamer2368 Mar 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 0 additions & 3 deletions include/linelast_solver.hpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should also be true for the Zero Neumann condition.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you sure about this? zero Neumann condition will be filtered out at the previous if clause, both previous and current version.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No I think that you're right. The examples all seem to work well now.

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
Open in desktop
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
Loading