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

Hanging nodes constraint fix for linear problems #48

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Hanging nodes constraint fix for linear problems #48

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
08fa5fa
Corrected hanging node constraints
bsumirak Jun 27, 2019
3f5d44d
Fix for bad implementation of hanging node constraints for linear pro…
bsumirak Jan 10, 2021
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
5 changes: 4 additions & 1 deletion ugbase/bridge/disc_bridges/algebra_bridge.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -123,6 +123,8 @@ static void Algebra(Registry& reg, string parentGroup)
.add_method("assemble_defect", static_cast<void (T::*)(vector_type&, const vector_type&, const GridLevel&)>(&T::assemble_defect),"", "d(u)#u#GridLevel", "Assembles Defect on grid level")
.add_method("assemble_linear", static_cast<void (T::*)(matrix_type&, vector_type&)>(&T::assemble_linear),"", "A#b", "Assembles Matrix and rhs on surface grid.")
.add_method("assemble_linear", static_cast<void (T::*)(matrix_type&, vector_type&, const GridLevel&)>(&T::assemble_linear),"", "A#b#GridLevel", "Assembles Matrix and rhs on grid level.")
.add_method("assemble_linear", static_cast<void (T::*)(matrix_type&, vector_type&, const vector_type&)>(&T::assemble_linear),"", "A#b#u", "Assembles Matrix and rhs on surface grid.")
.add_method("assemble_linear", static_cast<void (T::*)(matrix_type&, vector_type&, const vector_type&, const GridLevel&)>(&T::assemble_linear),"", "A#b#u#GridLevel", "Assembles Matrix and rhs on grid level.")
.add_method("assemble_rhs", static_cast<void (T::*)(vector_type&, const vector_type&)>(&T::assemble_rhs),"", "rhs#u", "assembles right-hand side on surface grid")
.add_method("assemble_rhs", static_cast<void (T::*)(vector_type&, const vector_type&, const GridLevel&)>(&T::assemble_rhs),"", "rhs#u#GridLevel", "assembles right-hand side on grid level")
.add_method("assemble_rhs", static_cast<void (T::*)(vector_type&)>(&T::assemble_rhs),"", "rhs#GridLevel", "assembles right-hand side on surface grid for linear case")
Expand Down Expand Up @@ -256,7 +258,8 @@ static void Algebra(Registry& reg, string parentGroup)
.add_method("set_discretization", &T::set_discretization)
.add_method("set_level", &T::set_level)
.add_method("set_dirichlet_values", &T::set_dirichlet_values)
.add_method("init_op_and_rhs", &T::init_op_and_rhs)
.add_method("init_op_and_rhs", static_cast<void (T::*) (vector_type&) >(&T::init_op_and_rhs))
.add_method("init_op_and_rhs", static_cast<void (T::*) (vector_type&, const vector_type&) >(&T::init_op_and_rhs))
.add_method("level", &T::level)
.set_construct_as_smart_pointer(true);
reg.add_class_to_group(name, "AssembledLinearOperator", tag);
Expand Down
1 change: 1 addition & 0 deletions ugbase/bridge/disc_bridges/domain_disc_bridge.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,7 @@ static void DomainAlgebra(Registry& reg, string grp)
.add_method("remove", static_cast<void (T::*)(SmartPtr<IElemDisc<TDomain> >)>(&T::remove), "", "Element Discretization")
.add_method("add", static_cast<void (T::*)(SmartPtr<IDiscretizationItem<TDomain, TAlgebra> >)>(&T::add), "", "DiscItem")
.add_method("assemble_linear", static_cast<void (T::*)(typename TAlgebra::matrix_type&, GridFunction<TDomain, TAlgebra>&)>(&T::assemble_linear))
.add_method("assemble_linear", static_cast<void (T::*)(typename TAlgebra::matrix_type&, GridFunction<TDomain, TAlgebra>&, const GridFunction<TDomain, TAlgebra>&)>(&T::assemble_linear))
.add_method("assemble_rhs", static_cast<void (T::*)(typename TAlgebra::vector_type&, GridFunction<TDomain, TAlgebra>&)>(&T::assemble_rhs))
.add_method("assemble_rhs", static_cast<void (T::*)(GridFunction<TDomain, TAlgebra>&)>(&T::assemble_rhs))
.add_method("adjust_solution", static_cast<void (T::*)(GridFunction<TDomain, TAlgebra>&)>(&T::adjust_solution))
Expand Down
2 changes: 1 addition & 1 deletion .../operator/preconditioner/projected_gauss_seidel/obstacles/obstacle_constraint_interface.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -200,7 +200,7 @@ class IObstacleConstraint:
};

/// sets unity rows in A and dirichlet values in right-hand side b
void adjust_linear(matrix_type& A, vector_type& b,
void adjust_linear(matrix_type& A, vector_type& b, const vector_type& u,
ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0){
UG_LOG("IObstacleConstraint::adjust_linear() \n");
};
Expand Down
14 changes: 12 additions & 2 deletions ugbase/lib_disc/assemble_interface.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -151,9 +151,19 @@ class IAssemble
* \param[out] b Right-Hand-Side
* \param[in] gl Grid Level
*/
virtual void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl) = 0;
virtual void assemble_linear(matrix_type& A, vector_type& b, const vector_type& u, const GridLevel& gl) = 0;
void assemble_linear(matrix_type& A, vector_type& b, const GridLevel& gl)
{
UG_LOGN("HINT: assemble_linear without given solution is deprecated.\n"
" Supply the solution as it will enter the solver as a third argument.");
// Use the rhs as solution. This will reproduce previous behavior
// when entries of b are assigned the corresponding entries of u.
assemble_linear(A, b, b, gl);
}
void assemble_linear(matrix_type& A, vector_type& b, const vector_type& u)
{assemble_linear(A, b, u, GridLevel());}
void assemble_linear(matrix_type& A, vector_type& b)
{assemble_linear(A,b, GridLevel());}
{assemble_linear(A, b, GridLevel());}

/// assembles rhs
virtual void assemble_rhs(vector_type& rhs, const vector_type& u, const GridLevel& gl) = 0;
Expand Down
3 changes: 3 additions & 0 deletions ugbase/lib_disc/operator/linear_operator/assembled_linear_operator.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,9 @@ class AssembledLinearOperator :
/// initializes the operator and assembles the passed rhs vector
void init_op_and_rhs(vector_type& b);

/// initializes the operator and assembles the passed rhs vector
void init_op_and_rhs(vector_type& b, const vector_type& u);

/// compute d = J(u)*c (here, J(u) is a Matrix)
virtual void apply(vector_type& d, const vector_type& c);

Expand Down
35 changes: 28 additions & 7 deletions ugbase/lib_disc/operator/linear_operator/assembled_linear_operator_impl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,9 @@ AssembledLinearOperator<TAlgebra>::init_op_and_rhs(vector_type& b)
{
// todo: check that assembling is linear

UG_LOGN("Hint: init_op_and_rhs without solution vector is deprecated.\n"
"Supply the solution vector (as it will enter the solver) as second argument.");

if(m_spAss.invalid())
UG_THROW("AssembledLinearOperator: Assembling routine not set.");

Expand All @@ -99,6 +102,24 @@ AssembledLinearOperator<TAlgebra>::init_op_and_rhs(vector_type& b)
" Cannot assemble Matrix and Rhs.");
}

// Initialize the operator
template <typename TAlgebra>
void
AssembledLinearOperator<TAlgebra>::init_op_and_rhs(vector_type& b, const vector_type& u)
{
// todo: check that assembling is linear

if(m_spAss.invalid())
UG_THROW("AssembledLinearOperator: Assembling routine not set.");

// assemble matrix and rhs in one loop
try{
m_spAss->assemble_linear(*this, b, u, m_gridLevel);
}
UG_CATCH_THROW("AssembledLinearOperator::init_op_and_rhs:"
" Cannot assemble Matrix and Rhs.");
}

template <typename TAlgebra>
void
AssembledLinearOperator<TAlgebra>::apply(vector_type& d, const vector_type& c)
Expand Down Expand Up @@ -170,20 +191,20 @@ void AssembleLinearOperatorRhsAndSolution
{
ASS_PROFILE_BEGIN(ASS_AssembleLinearOperatorRhsAndSolution);

// initialize operator
ASS_PROFILE_BEGIN(ASS_InitOperatorAndRhs);
try{
op.init_op_and_rhs(b);
}UG_CATCH_THROW("Cannot init the operator (assembling failed).");
ASS_PROFILE_END();

// sets the dirichlet values in the solution
ASS_PROFILE_BEGIN(ASS_SetDirValues);
try{
op.set_dirichlet_values(u);
} UG_CATCH_THROW("Cannot set the dirichlet values in the solution.");
ASS_PROFILE_END();

// initialize operator
ASS_PROFILE_BEGIN(ASS_InitOperatorAndRhs);
try{
op.init_op_and_rhs(b, u);
}UG_CATCH_THROW("Cannot init the operator (assembling failed).");
ASS_PROFILE_END();

ASS_PROFILE_END();
}

Expand Down
2 changes: 1 addition & 1 deletion ugbase/lib_disc/operator/linear_operator/nested_iteration/nested_iteration_impl.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -262,7 +262,7 @@ bool NestedIterationSolver<TDomain,TAlgebra>::apply(vector_type& u)

// Assemble.
NESTED_ITER_PROFILE_BEGIN(NestedIterationAssemble);
m_spAss->assemble_linear(*m_J, *spD, surfGridLevel); // todo: replace for non-linear problems
m_spAss->assemble_linear(*m_J, *spD, u, surfGridLevel); // todo: replace for non-linear problems
m_spAss->adjust_solution(u, surfGridLevel);
NESTED_ITER_PROFILE_END();

Expand Down
2 changes: 1 addition & 1 deletion ugbase/lib_disc/spatial_disc/constraints/constraint_interface.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ class IConstraint
const std::vector<number>* vScaleStiff = NULL) = 0;

/// adapts matrix and rhs (linear case) to enforce constraints
virtual void adjust_linear(matrix_type& mat, vector_type& rhs,
virtual void adjust_linear(matrix_type& mat, vector_type& rhs, const vector_type& u,
ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0) = 0;

/// adapts a rhs to enforce constraints
Expand Down
14 changes: 4 additions & 10 deletions ugbase/lib_disc/spatial_disc/constraints/continuity_constraints/p1_continuity_constraints.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ class SymP1Constraints
typedef typename algebra_type::vector_type vector_type;

public:
SymP1Constraints() : IDomainConstraint<TDomain, TAlgebra>(), m_bAssembleLinearProblem(false) {}
SymP1Constraints() : IDomainConstraint<TDomain, TAlgebra>() {}
virtual ~SymP1Constraints() {}

virtual int type() const {return CT_HANGING;}
Expand All @@ -80,7 +80,7 @@ class SymP1Constraints
void adjust_rhs(vector_type& rhs, const vector_type& u,
ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0);

void adjust_linear(matrix_type& mat, vector_type& rhs,
void adjust_linear(matrix_type& mat, vector_type& rhs, const vector_type& u,
ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0);

void adjust_solution(vector_type& u, ConstSmartPtr<DoFDistribution> dd,
Expand All @@ -104,9 +104,6 @@ class SymP1Constraints
int type,
number time = 0.0
);

protected:
bool m_bAssembleLinearProblem;
};


Expand All @@ -129,7 +126,7 @@ class OneSideP1Constraints
typedef IDomainConstraint<TDomain, TAlgebra> base_type;

public:
OneSideP1Constraints() : IDomainConstraint<TDomain, TAlgebra>(), m_bAssembleLinearProblem(false) {}
OneSideP1Constraints() : IDomainConstraint<TDomain, TAlgebra>() {}
virtual ~OneSideP1Constraints() {}

virtual int type() const {return CT_HANGING;}
Expand All @@ -148,7 +145,7 @@ class OneSideP1Constraints
void adjust_rhs(vector_type& rhs, const vector_type& u,
ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0);

void adjust_linear(matrix_type& mat, vector_type& rhs,
void adjust_linear(matrix_type& mat, vector_type& rhs, const vector_type& u,
ConstSmartPtr<DoFDistribution> dd, int type, number time = 0.0);

void adjust_solution(vector_type& u, ConstSmartPtr<DoFDistribution> dd,
Expand All @@ -172,9 +169,6 @@ class OneSideP1Constraints
int type,
number time = 0.0
);

protected:
bool m_bAssembleLinearProblem;
};

}; // namespace ug
Expand Down
Loading