NormsOperators.cpp

Go to the documentation of this file.

/** \file NormsOperators.cpp
 
\brief User data operators for calculating norms and differences between
fields
 
*/
 
namespace MoFEM {
 
OpCalcNormL2Tensor0::OpCalcNormL2Tensor0(
    boost::shared_ptr<VectorDouble> data_ptr, SmartPetscObj<Vec> data_vec,
    const int index, boost::shared_ptr<VectorDouble> diff_data_ptr,
    boost::shared_ptr<Range> ent_ptr)
    : ForcesAndSourcesCore::UserDataOperator(
          NOSPACE, ForcesAndSourcesCore::UserDataOperator::OPSPACE),
      dataPtr(data_ptr), dataVec(data_vec), iNdex(index),
      diffDataPtr(diff_data_ptr), entsPtr(ent_ptr) {
  if (!dataPtr)
    THROW_MESSAGE("Pointer is not set");
  if (!diffDataPtr)
    diffDataPtr = dataPtr;
}
 
MoFEMErrorCode OpCalcNormL2Tensor0::doWork(int side, EntityType type,
                                           EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
 
  // check if entity is in the range
  if (entsPtr) {
    if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
      MoFEMFunctionReturnHot(0);
  }
 
  // calculate the difference between data pointers and save them in diffDataPtr
  if (dataPtr != diffDataPtr)
    *diffDataPtr -= *dataPtr;
 
  // get number of integration points
  const auto nb_integration_points = getGaussPts().size2();
  // get element volume
  const double vol = getMeasure();
  // get integration weights
  auto t_w = getFTensor0IntegrationWeight();
  // get values
  auto t_data = getFTensor0FromVec(*diffDataPtr);
  // initialise double to store norm values
  double norm_on_element = 0.;
  // loop over integration points
  for (int gg = 0; gg != nb_integration_points; gg++) {
    // add to element norm
    norm_on_element += t_w * t_data * t_data;
    // move to another integration weight
    ++t_w;
    // move to another data values
    ++t_data;
  }
  // scale with volume of the element
  norm_on_element *= vol;
  // add to dataVec at iNdex position
  CHKERR VecSetValue(dataVec, iNdex, norm_on_element, ADD_VALUES);
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM>
OpCalcNormL2Tensor1<DIM>::OpCalcNormL2Tensor1(
    boost::shared_ptr<MatrixDouble> data_ptr, SmartPetscObj<Vec> data_vec,
    const int index, boost::shared_ptr<MatrixDouble> diff_data_ptr,
    boost::shared_ptr<Range> ent_ptr)
    : ForcesAndSourcesCore::UserDataOperator(
          NOSPACE, ForcesAndSourcesCore::UserDataOperator::OPSPACE),
      dataPtr(data_ptr), dataVec(data_vec), iNdex(index),
      diffDataPtr(diff_data_ptr), entsPtr(ent_ptr) {
  if (!dataPtr)
    THROW_MESSAGE("Pointer is not set");
  if (!diffDataPtr)
    diffDataPtr = dataPtr;
}
 
template <int DIM>
MoFEMErrorCode
OpCalcNormL2Tensor1<DIM>::doWork(int side, EntityType type,
                                 EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
 
  // check if entity is in the range
  if (entsPtr) {
    if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
      MoFEMFunctionReturnHot(0);
  }
 
  // calculate the difference between data pointers and save them in diffDataPtr
  if (dataPtr != diffDataPtr)
    *diffDataPtr -= *dataPtr;
 
  // Declare FTensor index
  FTensor::Index<'i', DIM> i;
  // get number of integration points
  const auto nb_integration_points = getGaussPts().size2();
  // get element volume
  const double vol = getMeasure();
  // get integration weights
  auto t_w = getFTensor0IntegrationWeight();
  // get vector values
#ifndef NDEBUG
  if (diffDataPtr->size2() != DIM)
    THROW_MESSAGE("OpCalcNormL2Tensor1: wrong number of matrix columns");
#endif
  auto t_data = getFTensor1FromPtr<DIM, DIM>(diffDataPtr->data().data());
  // initialise double to store norm values
  double norm_on_element = 0.;
  // loop over integration points
  for (int gg = 0; gg != nb_integration_points; ++gg) {
    // add to element norm
    norm_on_element += t_w * (t_data(i) * t_data(i));
    // move to another integration weight
    ++t_w;
    // move to another data values
    ++t_data;
  }
  // scale with volume of the element
  norm_on_element *= vol;
  // add to dataVec at iNdex position
  CHKERR VecSetValue(dataVec, iNdex, norm_on_element, ADD_VALUES);
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM_1, int DIM_2>
OpCalcNormL2Tensor2<DIM_1, DIM_2>::OpCalcNormL2Tensor2(
    boost::shared_ptr<MatrixDouble> data_ptr, SmartPetscObj<Vec> data_vec,
    const int index, boost::shared_ptr<MatrixDouble> diff_data_ptr,
    boost::shared_ptr<Range> ent_ptr)
    : ForcesAndSourcesCore::UserDataOperator(
          NOSPACE, ForcesAndSourcesCore::UserDataOperator::OPSPACE),
      dataPtr(data_ptr), dataVec(data_vec), iNdex(index),
      diffDataPtr(diff_data_ptr), entsPtr(ent_ptr) {
  if (!dataPtr)
    THROW_MESSAGE("Pointer is not set");
  if (!diffDataPtr)
    diffDataPtr = dataPtr;
}
 
template <int DIM_1, int DIM_2>
MoFEMErrorCode
OpCalcNormL2Tensor2<DIM_1, DIM_2>::doWork(int side, EntityType type,
                                          EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
 
  // check if entity is in the range
  if (entsPtr) {
    if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
      MoFEMFunctionReturnHot(0);
  }
 
  // calculate the difference between data pointers and save them in diffDataPtr
  if (dataPtr != diffDataPtr)
    *diffDataPtr -= *dataPtr;
 
  // Declare FTensor index
  FTensor::Index<'i', DIM_1> i;
  FTensor::Index<'j', DIM_2> j;
  // get number of integration points
  const auto nb_integration_points = getGaussPts().size2();
  // get element volume
  const double vol = getMeasure();
  // get integration weights
  auto t_w = getFTensor0IntegrationWeight();
  // get vector values
  auto t_data = getFTensor2FromMat<DIM_1, DIM_2>(*diffDataPtr);
  // initialise double to store norm values
  double norm_on_element = 0.;
  // loop over integration points
  for (int gg = 0; gg != nb_integration_points; gg++) {
    // add to element norm
    norm_on_element += t_w * (t_data(i, j) * t_data(i, j));
    // move to another integration weight
    ++t_w;
    // move to another data values
    ++t_data;
  }
  // scale with volume of the element
  norm_on_element *= vol;
  // add to dataVec at iNdex position
  CHKERR VecSetValue(dataVec, iNdex, norm_on_element, ADD_VALUES);
 
  MoFEMFunctionReturn(0);
}
 
template <int DIM>
MoFEMErrorCode
OpCalcNormL2Tensor2Symmetric<DIM>::doWork(int side, EntityType type,
                                          EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
 
  // check if entity is in the range
  if (this->entsPtr) {
    if (this->entsPtr->find(this->getFEEntityHandle()) == this->entsPtr->end())
      MoFEMFunctionReturnHot(0);
  }
 
  // calculate the difference between data pointers and save them in diffDataPtr
  if (this->dataPtr != this->diffDataPtr)
    *this->diffDataPtr -= *this->dataPtr;
 
  // Declare FTensor index
  FTensor::Index<'i', DIM> i;
  FTensor::Index<'j', DIM> j;
  // get number of integration points
  const auto nb_integration_points = this->getGaussPts().size2();
  // get element volume
  const double vol = this->getMeasure();
  // get integration weights
  auto t_w = this->getFTensor0IntegrationWeight();
  // get symmetric tensor values
  auto t_data = getFTensor2SymmetricFromMat<DIM>(*this->diffDataPtr);
  // initialise double to store norm values
  double norm_on_element = 0.;
  // loop over integration points
  for (int gg = 0; gg != nb_integration_points; gg++) {
    // add Frobenius norm contribution, including both off-diagonal entries
    norm_on_element += t_w * (t_data(i, j) * t_data(i, j));
    // move to another integration weight
    ++t_w;
    // move to another data values
    ++t_data;
  }
  // scale with volume of the element
  norm_on_element *= vol;
  // add to dataVec at iNdex position
  CHKERR VecSetValue(this->dataVec, this->iNdex, norm_on_element, ADD_VALUES);
 
  MoFEMFunctionReturn(0);
}
 
OpGetTensor0fromFunc::OpGetTensor0fromFunc(
    boost::shared_ptr<VectorDouble> data_ptr, ScalarFun scalar_function)
    : ForcesAndSourcesCore::UserDataOperator(
          NOSPACE, ForcesAndSourcesCore::UserDataOperator::OPSPACE),
      dataPtr(data_ptr), sFunc(scalar_function) {}
 
MoFEMErrorCode OpGetTensor0fromFunc::doWork(int side, EntityType type,
                                            EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
 
  // get number of integration points
  const int nb_integration_pts = getGaussPts().size2();
  // resize dataPtr to store values at integration points
  dataPtr->resize(nb_integration_pts);
  // get values
  auto t_val = getFTensor0FromVec(*dataPtr);
 
  // get coordinates at integration point
  auto t_coords = getFTensor1CoordsAtGaussPts();
  // loop over integration points
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    // set value at integration point to value of scalar function at the
    // coordinates
    t_val = sFunc(t_coords(0), t_coords(1), t_coords(2));
    // move to another value
    ++t_val;
    // move to another coordinates
    ++t_coords;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int TENS_DIM, int FUNC_DIM>
OpGetTensor1fromFunc<TENS_DIM, FUNC_DIM>::OpGetTensor1fromFunc(
    boost::shared_ptr<MatrixDouble> data_ptr, VectorFunc vector_function)
    : ForcesAndSourcesCore::UserDataOperator(
          NOSPACE, ForcesAndSourcesCore::UserDataOperator::OPSPACE),
      dataPtr(data_ptr), vFunc(vector_function) {}
 
template <int TENS_DIM, int FUNC_DIM>
MoFEMErrorCode OpGetTensor1fromFunc<TENS_DIM, FUNC_DIM>::doWork(
    int side, EntityType type, EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
  if (TENS_DIM > FUNC_DIM) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Output tensor dimension TENS_DIM (%d) is larger than FUNC_DIM "
            "(%d)",
            TENS_DIM, FUNC_DIM);
  }
  // Copy only the components requested by the output tensor. The function
  // buffer can be larger, e.g. a 2D tensor read from a 3D function buffer.
  FTensor::Index<'i', TENS_DIM> i;
  // get number of integration points
  const int nb_integration_pts = getGaussPts().size2();
  // resize dataPtr to store values at integration points
  using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
  auto get_values_at_pts =
      MatrixSizeHelper<GetFTensor1FromMatType<TENS_DIM, -1, DL>, DL>::size(
          *dataPtr, nb_integration_pts);
  dataPtr->clear();
  // get values
  auto t_val = get_values_at_pts();
 
  // get coordinates at integration point
  auto t_coords = getFTensor1CoordsAtGaussPts();
  // loop over integration points
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    // get function values of vector function at the coordinates
    auto func_val = vFunc(t_coords(0), t_coords(1), t_coords(2));
    // translate function values to tensor
    auto t_func_val = getFTensor1FromArray<FUNC_DIM, FUNC_DIM>(func_val);
    // set values at integration point to function values at the coordinates
    t_val(i) = t_func_val(i);
 
    // move to another value
    ++t_val;
    // move to another coordinates
    ++t_coords;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int TENS_DIM, int FUNC_DIM>
OpGetTensor2fromFunc<TENS_DIM, FUNC_DIM>::OpGetTensor2fromFunc(
    boost::shared_ptr<MatrixDouble> data_ptr, MatrixFunc matrix_function)
    : ForcesAndSourcesCore::UserDataOperator(
          NOSPACE, ForcesAndSourcesCore::UserDataOperator::OPSPACE),
      dataPtr(data_ptr), mFunc(matrix_function) {}
 
template <int TENS_DIM, int FUNC_DIM>
MoFEMErrorCode OpGetTensor2fromFunc<TENS_DIM, FUNC_DIM>::doWork(
    int side, EntityType type, EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
  if (TENS_DIM > FUNC_DIM) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Output tensor dimension TENS_DIM (%d) is larger than FUNC_DIM "
            "(%d)",
            TENS_DIM, FUNC_DIM);
  }
  // Copy only the components requested by the output tensor. The function
  // buffer can be larger, e.g. a 2D tensor read from a 3D function buffer.
  FTensor::Index<'i', TENS_DIM> i;
  FTensor::Index<'j', TENS_DIM> j;
  // get number of integration points
  const int nb_integration_pts = getGaussPts().size2();
  // resize dataPtr to store values at integration points
  dataPtr->resize(TENS_DIM * TENS_DIM, nb_integration_pts);
  // get values
  auto t_val = getFTensor2FromMat<TENS_DIM, TENS_DIM>(*(dataPtr));
 
  // get coordinates at integration point
  auto t_coords = getFTensor1CoordsAtGaussPts();
  // loop over integration points
  for (int gg = 0; gg != nb_integration_pts; ++gg) {
    // get function values of matrix function at the coordinates
    auto func_val = mFunc(t_coords(0), t_coords(1), t_coords(2));
    // translate function values to tensor
    auto t_func_val =
        getFTensor2FromArray<FUNC_DIM, FUNC_DIM, FUNC_DIM>(func_val, 0);
    // set values at integration point to function values at the coordinates
    t_val(i, j) = t_func_val(i, j);
 
    // move to another value
    ++t_val;
    // move to another coordinates
    ++t_coords;
  }
 
  MoFEMFunctionReturn(0);
}
 
template <int TENS_DIM, int FUNC_DIM>
MoFEMErrorCode OpGetTensor2SymmetricfromFunc<TENS_DIM, FUNC_DIM>::doWork(
    int side, EntityType type, EntitiesFieldData::EntData &data) {
  MoFEMFunctionBegin;
 
  if constexpr (TENS_DIM > FUNC_DIM) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Output tensor dimension TENS_DIM (%d) is larger than FUNC_DIM "
            "(%d)",
            TENS_DIM, FUNC_DIM);
  }
 
  const auto nb_gauss_pts = this->getGaussPts().size2();
  using DL = DataLayoutTraits<DataLayout::GaussByCoeffs>;
  auto get_values_at_gauss_pts =
      MatrixSizeHelper<GetFTensor2SymmetricFromMatType<TENS_DIM, -1, DL>,
                       DL>::size(*this->dataPtr, nb_gauss_pts);
  this->dataPtr->clear();
 
  auto t_values_at_gauss_pts = get_values_at_gauss_pts();
  auto t_coords = this->getFTensor1CoordsAtGaussPts();
  FTensor::Index<'i', TENS_DIM> i;
  FTensor::Index<'j', TENS_DIM> j;
 
  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
    auto values = this->mFunc(t_coords(0), t_coords(1), t_coords(2));
    auto t_values = getFTensor2SymmetricFromMat<
        FUNC_DIM, -1, DataLayoutTraits<DataLayout::CoeffsByGauss>>(values);
    t_values_at_gauss_pts(i, j) = t_values(i, j);
    ++t_values_at_gauss_pts;
    ++t_coords;
  }
 
  MoFEMFunctionReturn(0);
}
 
template struct OpCalcNormL2Tensor1<2>;
template struct OpCalcNormL2Tensor1<3>;
template struct OpCalcNormL2Tensor2<2, 2>;
template struct OpCalcNormL2Tensor2<3, 3>;
template struct OpCalcNormL2Tensor2Symmetric<2>;
template struct OpCalcNormL2Tensor2Symmetric<3>;
template struct OpGetTensor1fromFunc<2, 2>;
template struct OpGetTensor1fromFunc<2, 3>;
template struct OpGetTensor1fromFunc<3, 3>;
template struct OpGetTensor2fromFunc<2, 2>;
template struct OpGetTensor2fromFunc<2, 3>;
template struct OpGetTensor2fromFunc<3, 3>;
template struct OpGetTensor2SymmetricfromFunc<2, 2>;
template struct OpGetTensor2SymmetricfromFunc<2, 3>;
template struct OpGetTensor2SymmetricfromFunc<3, 3>;
 
} // namespace MoFEM