<commit_before>//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "InitialConditionTempl.h"
#include "FEProblem.h"
#include "Assembly.h"
#include "MooseVariableFE.h"
#include "SystemBase.h"

#include "libmesh/fe_interface.h"
#include "libmesh/quadrature.h"

template <>
InputParameters
validParams<InitialConditionTempl<Real>>()
{
  return InitialConditionTempl<Real>::validParams();
}
template <>
InputParameters
validParams<InitialConditionTempl<RealVectorValue>>()
{
  return InitialConditionTempl<RealVectorValue>::validParams();
}

template <>
InputParameters
validParams<InitialConditionTempl<RealEigenVector>>()
{
  return InitialConditionTempl<RealEigenVector>::validParams();
}

template <typename T>
InputParameters
InitialConditionTempl<T>::validParams()
{
  InputParameters params = InitialConditionBase::validParams();
  return params;
}

template <typename T>
InitialConditionTempl<T>::InitialConditionTempl(const InputParameters & parameters)
  : InitialConditionBase(parameters),
    _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
    _tid(getParam<THREAD_ID>("_tid")),
    _assembly(_fe_problem.assembly(_tid)),
    _t(_fe_problem.time()),
    _coord_sys(_assembly.coordSystem()),
    _var(_sys.getFieldVariable<T>(parameters.get<THREAD_ID>("_tid"),
                                  parameters.get<VariableName>("variable"))),
    _current_elem(_var.currentElem()),
    _current_node(nullptr),
    _qp(0),
    _fe_type(_var.feType())
{
}

template <typename T>
InitialConditionTempl<T>::~InitialConditionTempl()
{
}

template <typename T>
void
InitialConditionTempl<T>::compute()
{
  // -- NOTE ----
  // The following code is a copy from libMesh project_vector.C plus it adds some features, so we
  // can couple variable values
  // and we also do not call any callbacks, but we use our initial condition system directly.
  // ------------

  // The dimension of the current element
  _dim = _current_elem->dim();
  // The element type
  const ElemType elem_type = _current_elem->type();
  // The number of nodes on the new element
  const unsigned int n_nodes = _current_elem->n_nodes();

  // Get FE objects of the appropriate type
  // We cannot use the FE object in Assembly, since the following code is messing with the
  // quadrature rules
  // for projections and would screw it up. However, if we implement projections from one mesh to
  // another,
  // this code should use that implementation.
  std::unique_ptr<FEBaseType> fe(FEBaseType::build(_dim, _fe_type));

  // Prepare variables for projection
  std::unique_ptr<QBase> qrule(_fe_type.default_quadrature_rule(_dim));
  std::unique_ptr<QBase> qedgerule(_fe_type.default_quadrature_rule(1));
  std::unique_ptr<QBase> qsiderule(_fe_type.default_quadrature_rule(_dim - 1));

  // The values of the shape functions at the quadrature points
  _phi = &fe->get_phi();

  // The gradients of the shape functions at the quadrature points on the child element.
  _dphi = nullptr;

  _cont = fe->get_continuity();

  if (_cont == C_ONE)
  {
    const std::vector<std::vector<GradientShapeType>> & ref_dphi = fe->get_dphi();
    _dphi = &ref_dphi;
  }

  // The Jacobian * quadrature weight at the quadrature points
  _JxW = &fe->get_JxW();
  // The XYZ locations of the quadrature points
  _xyz_values = &fe->get_xyz();

  // Update the DOF indices for this element based on the current mesh
  _var.prepareIC();
  _dof_indices = _var.dofIndices();

  // The number of DOFs on the element
  const unsigned int n_dofs = _dof_indices.size();
  if (n_dofs == 0)
    return;

  // Fixed vs. free DoFs on edge/face projections
  _dof_is_fixed.clear();
  _dof_is_fixed.resize(n_dofs, false);
  _free_dof.clear();
  _free_dof.resize(n_dofs, 0);

  // Zero the interpolated values
  _Ue.resize(n_dofs);
  _Ue.zero();

  DenseVector<char> mask(n_dofs, true);

  // In general, we need a series of
  // projections to ensure a unique and continuous
  // solution.  We start by interpolating nodes, then
  // hold those fixed and project edges, then
  // hold those fixed and project faces, then
  // hold those fixed and project interiors

  // Interpolate node values first
  _current_dof = 0;

  for (_n = 0; _n != n_nodes; ++_n)
  {
    _nc = FEInterface::n_dofs_at_node(_dim, _fe_type, elem_type, _n);

    // for nodes that are in more than one subdomain, only compute the initial
    // condition once on the lowest numbered block
    auto curr_node = _current_elem->node_ptr(_n);
    const auto & block_ids = _sys.mesh().getNodeBlockIds(*curr_node);
    auto priority_block = *block_ids.begin();

    if (!hasBlocks(priority_block) && _var.isNodal())
    {
      for (decltype(_nc) i = 0; i < _nc; ++i)
      {
        mask(_current_dof) = false;
        _current_dof++;
      }
      continue;
    }

    // FIXME: this should go through the DofMap,
    // not duplicate _dof_indices code badly!
    if (!_current_elem->is_vertex(_n))
    {
      _current_dof += _nc;
      continue;<commit_msg>fix priority block selection<commit_after>//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "InitialConditionTempl.h"
#include "FEProblem.h"
#include "Assembly.h"
#include "MooseVariableFE.h"
#include "SystemBase.h"

#include "libmesh/fe_interface.h"
#include "libmesh/quadrature.h"

template <>
InputParameters
validParams<InitialConditionTempl<Real>>()
{
  return InitialConditionTempl<Real>::validParams();
}
template <>
InputParameters
validParams<InitialConditionTempl<RealVectorValue>>()
{
  return InitialConditionTempl<RealVectorValue>::validParams();
}

template <>
InputParameters
validParams<InitialConditionTempl<RealEigenVector>>()
{
  return InitialConditionTempl<RealEigenVector>::validParams();
}

template <typename T>
InputParameters
InitialConditionTempl<T>::validParams()
{
  InputParameters params = InitialConditionBase::validParams();
  return params;
}

template <typename T>
InitialConditionTempl<T>::InitialConditionTempl(const InputParameters & parameters)
  : InitialConditionBase(parameters),
    _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
    _tid(getParam<THREAD_ID>("_tid")),
    _assembly(_fe_problem.assembly(_tid)),
    _t(_fe_problem.time()),
    _coord_sys(_assembly.coordSystem()),
    _var(_sys.getFieldVariable<T>(parameters.get<THREAD_ID>("_tid"),
                                  parameters.get<VariableName>("variable"))),
    _current_elem(_var.currentElem()),
    _current_node(nullptr),
    _qp(0),
    _fe_type(_var.feType())
{
}

template <typename T>
InitialConditionTempl<T>::~InitialConditionTempl()
{
}

template <typename T>
void
InitialConditionTempl<T>::compute()
{
  // -- NOTE ----
  // The following code is a copy from libMesh project_vector.C plus it adds some features, so we
  // can couple variable values
  // and we also do not call any callbacks, but we use our initial condition system directly.
  // ------------

  // The dimension of the current element
  _dim = _current_elem->dim();
  // The element type
  const ElemType elem_type = _current_elem->type();
  // The number of nodes on the new element
  const unsigned int n_nodes = _current_elem->n_nodes();

  // Get FE objects of the appropriate type
  // We cannot use the FE object in Assembly, since the following code is messing with the
  // quadrature rules
  // for projections and would screw it up. However, if we implement projections from one mesh to
  // another,
  // this code should use that implementation.
  std::unique_ptr<FEBaseType> fe(FEBaseType::build(_dim, _fe_type));

  // Prepare variables for projection
  std::unique_ptr<QBase> qrule(_fe_type.default_quadrature_rule(_dim));
  std::unique_ptr<QBase> qedgerule(_fe_type.default_quadrature_rule(1));
  std::unique_ptr<QBase> qsiderule(_fe_type.default_quadrature_rule(_dim - 1));

  // The values of the shape functions at the quadrature points
  _phi = &fe->get_phi();

  // The gradients of the shape functions at the quadrature points on the child element.
  _dphi = nullptr;

  _cont = fe->get_continuity();

  if (_cont == C_ONE)
  {
    const std::vector<std::vector<GradientShapeType>> & ref_dphi = fe->get_dphi();
    _dphi = &ref_dphi;
  }

  // The Jacobian * quadrature weight at the quadrature points
  _JxW = &fe->get_JxW();
  // The XYZ locations of the quadrature points
  _xyz_values = &fe->get_xyz();

  // Update the DOF indices for this element based on the current mesh
  _var.prepareIC();
  _dof_indices = _var.dofIndices();

  // The number of DOFs on the element
  const unsigned int n_dofs = _dof_indices.size();
  if (n_dofs == 0)
    return;

  // Fixed vs. free DoFs on edge/face projections
  _dof_is_fixed.clear();
  _dof_is_fixed.resize(n_dofs, false);
  _free_dof.clear();
  _free_dof.resize(n_dofs, 0);

  // Zero the interpolated values
  _Ue.resize(n_dofs);
  _Ue.zero();

  DenseVector<char> mask(n_dofs, true);

  // In general, we need a series of
  // projections to ensure a unique and continuous
  // solution.  We start by interpolating nodes, then
  // hold those fixed and project edges, then
  // hold those fixed and project faces, then
  // hold those fixed and project interiors

  // Interpolate node values first
  _current_dof = 0;

  for (_n = 0; _n != n_nodes; ++_n)
  {
    _nc = FEInterface::n_dofs_at_node(_dim, _fe_type, elem_type, _n);

    // for nodes that are in more than one subdomain, only compute the initial
    // condition once on the lowest numbered block
    auto curr_node = _current_elem->node_ptr(_n);
    const auto & block_ids = _sys.mesh().getNodeBlockIds(*curr_node);

    auto priority_block = *(block_ids.begin());
    for (auto id : block_ids)
      if (_var.hasBlocks(id))
      {
        priority_block = id;
        break;
      }

    if (!hasBlocks(priority_block) && _var.isNodal())
    {
      for (decltype(_nc) i = 0; i < _nc; ++i)
      {
        mask(_current_dof) = false;
        _current_dof++;
      }
      continue;
    }

    // FIXME: this should go through the DofMap,
    // not duplicate _dof_indices code badly!
    if (!_current_elem->is_vertex(_n))
    {
      _current_dof += _nc;
      continue;
{'fraction_non_alphanumeric': 0.11842351619546901, 'fraction_numerical': 0.0022467702677401234, 'mean_word_length': 3.8317503392130257, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 87, 'https://': 6, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3857602', 'n_tokens_mistral': 3465, 'n_tokens_neox': 3372, 'n_words': 919}