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.

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

First version of constraint boundary conditions #2145

Merged
merged 15 commits into from
Jun 27, 2018

Conversation

TomFischer
Copy link
Member

For discussion:

  • Output information via DBUG or INFO?
  • Open question: How to specify which getFlux function (for instance heat, mass or Darcy flux) should be used for the constraint calculation?

ProcessLib::createLocalAssemblers<
ConstraintDirichletBoundaryConditionLocalAssembler>(
_bulk_mesh.getDimension(), _bc_mesh.getElements(),
*_dof_table_boundary, 1, _local_assemblers, false, _integration_order,
Copy link
Member

Choose a reason for hiding this comment

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

How about to introduce two constants as:

    const int component = 1;
    const bool lower = false; 

to replace the numbers, 1. false, in the argument list?

Copy link
Member Author

Choose a reason for hiding this comment

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

Replacing the number makes the code more readable. Are you sure the suggested names are correct? What do think about the following replacements:

const int shape_function_order = 1;
const bool is_axisymmetric = false;

?

Copy link
Member

Choose a reason for hiding this comment

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

Sorry I've checked a wrong argument list. The names of shape_function_order and is_axisymmetric are correct.

unsigned const number_nodes = boundary_element->getNumberOfNodes();
for (unsigned i = 0; i < number_nodes; ++i)
{
unsigned const id = boundary_element->getNode(i)->getID();
Copy link
Member

Choose a reason for hiding this comment

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

getID() returns std::size_t.

Copy link
Member Author

Choose a reason for hiding this comment

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

Changed to auto.

class ConstraintDirichletBoundaryCondition final : public BoundaryCondition
{
public:
/// @param parameter used for setting the values for the boundary condition
Copy link
Member

Choose a reason for hiding this comment

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

It is better to capitalize the first letter of argument comment.

Copy link
Member Author

Choose a reason for hiding this comment

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

Capitalized.

}

// make the pairs (node id, value) unique
std::sort(tmp_bc_values.begin(), tmp_bc_values.end(),
Copy link
Member

Choose a reason for hiding this comment

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

The comment is about making vector unique, but here it's only sorted... ???

/// @param integration_order The order of the integration.
ConstraintDirichletBoundaryConditionLocalAssembler(
MeshLib::Element const& surface_element,
std::size_t /*const local_matrix_size*/,
Copy link
Member

Choose a reason for hiding this comment

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

Commenting 'const' changes function's signature...


IntegrationMethod const _integration_method;
std::size_t _bulk_element_id;
unsigned _bulk_face_id;
Copy link
Member

Choose a reason for hiding this comment

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

these _bulk_*_id are not const?

});

const int shape_function_order = 1;
const bool is_axisymmetric = false;
Copy link
Member

Choose a reason for hiding this comment

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

isn't this a property of the bc_mesh?

// check if the boundary element is active
if (_lower)
{
if (_flux_values[boundary_element->getID()] < _constraint_threshold)
Copy link
Member

Choose a reason for hiding this comment

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

A comparator before the for-loop would be better:

auto is_flux = [&](std::size_t const element_id) {
    return _lower ? _flux_values[element_id] < _constraint_threshold
       : _flux_values[element_id] > _constraint_threshold; }
for ( ... boundary elements ...)
{
   if (is_flux(boundary_element->getID())
   {
       continue;
   }
   ...
}

//! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__ConstraintDirichletBoundaryCondition__constraining_process_variable}
config.getConfigParameter<std::string>("constraining_process_variable");

const int process_id = 0; // assume monolithic implementation
Copy link
Member

Choose a reason for hiding this comment

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

Maybe this can be checked if the constraining_process is monolithic?

@TomFischer TomFischer force-pushed the ConstraintBCs branch 3 times, most recently from d71e80d to 7db4f7b Compare June 22, 2018 09:40
private:
MeshLib::Element const& _surface_element;

std::vector<double> _detJ_times_integralMeasure;
Copy link
Member

Choose a reason for hiding this comment

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

integration point weight could also be added (multiplied) in the same scalar...

// for correct results it is necessary to multiply the normal with -1
surface_element_normal *= -1;

std::size_t const n_integration_points =
Copy link
Member

Choose a reason for hiding this comment

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

I'd recommend to use 'auto' here... The return type might change, and is currently an 'unsigned int'.

ConstViscosityThermalConvection_pcs_0_ts_3_t_0.001010.vtu ConstViscosityThermalConvection_pcs_0_ts_3_t_0.001010.vtu p p 1e-15 1e-14
ConstViscosityThermalConvection_pcs_0_ts_4_t_0.101010.vtu ConstViscosityThermalConvection_pcs_0_ts_4_t_0.101010.vtu p p 1e-15 1e-14
ConstViscosityThermalConvection_pcs_0_ts_5_t_1.101010.vtu ConstViscosityThermalConvection_pcs_0_ts_5_t_1.101010.vtu p p 1e-15 1e-14
ConstViscosityThermalConvection_pcs_0_ts_6_t_10.000000.vtu ConstViscosityThermalConvection_pcs_0_ts_6_t_10.000000.vtu p p 1e-15 1e-14
Copy link
Member

Choose a reason for hiding this comment

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


// make the pairs (node id, value) unique
// first: sort the (node id, value) pairs according to the node id
std::sort(tmp_bc_values.begin(), tmp_bc_values.end(),
Copy link
Member

Choose a reason for hiding this comment

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

Sorry, but I still don't see where the tmp_bc_values are "made unique". There is a copy, but not a std::unique followed by erase....

Copy link
Member Author

Choose a reason for hiding this comment

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

I changed the comment.

GlobalVector const&)>
getFlux);

~ConstraintDirichletBoundaryCondition() = default;
Copy link
Member

Choose a reason for hiding this comment

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

Can be dropped, isn't it?

// At the moment (2018-04-26) the surface normal is not oriented
// according to the right hand rule
// for correct results it is necessary to multiply the normal with -1
surface_element_normal *= -1;
Copy link
Member

Choose a reason for hiding this comment

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

Computing a normalized surface normal for element/face could be a free function in MeshLib.

Eigen::Map<Eigen::RowVectorXd const>(bulk_flux.data(),
bulk_flux.size())
.dot(Eigen::Map<Eigen::RowVectorXd const>(
surface_element_normal.getCoords(), 3)));
Copy link
Member

Choose a reason for hiding this comment

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

Since the surface element normal does not change, it could be stored in the integration point data...

Copy link
Member Author

Choose a reason for hiding this comment

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

The normal doesn't change for all integration points, am I right? For this reason, I made the normal an attribute of the object and did not move it to the integration point data.

Copy link
Member

Choose a reason for hiding this comment

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

The current implementation is returning a constant normal vector for each element. This is not correct...

@endJunction endJunction merged commit e54815c into ufz:master Jun 27, 2018
@TomFischer TomFischer deleted the ConstraintBCs branch June 28, 2018 04:06
bilke pushed a commit to bilke/ogs that referenced this pull request Jul 2, 2018
First version of constraint boundary conditions
@chleh chleh mentioned this pull request Aug 13, 2018
15 tasks
@ogsbot
Copy link
Member

ogsbot commented Jun 19, 2020

OpenGeoSys development has been moved to GitLab.

See this pull request on GitLab.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
4 participants