v0.13.0
NavierStokesElement.cpp

Implementation of operators for fluid flowImplementation of operators for simulation of the fluid flow governed by Stokes and Navier-Stokes equations, and computation of the drag force on a fluid-solid inteface

/**
* \file NavierStokesElement.cpp
* \example NavierStokesElement.cpp
*
* \brief Implementation of operators for fluid flow
*
* Implementation of operators for simulation of the fluid flow governed by
* Stokes and Navier-Stokes equations, and computation of the drag
* force on a fluid-solid inteface
*/
/* This file is part of MoFEM.
* MoFEM is free software: you can redistribute it and/or modify it under
* the terms of the GNU Lesser General Public License as published by the
* Free Software Foundation, either version 3 of the License, or (at your
* option) any later version.
*
* MoFEM is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
* License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
#include <MoFEM.hpp>
using namespace MoFEM;
using namespace boost::numeric;
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_rhs_ptr,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs_ptr,
const std::string velocity_field, const std::string pressure_field,
boost::shared_ptr<CommonData> common_data, const EntityType type) {
for (auto &sit : common_data->setOfBlocksData) {
if (type == MBPRISM) {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
fe_lhs_ptr->getOpPtrVector().push_back(
new OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms());
fe_lhs_ptr->getOpPtrVector().push_back(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms());
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
}
fe_lhs_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
velocity_field, common_data->velPtr));
fe_lhs_ptr->getOpPtrVector().push_back(
common_data->gradVelPtr));
fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsDiagLin(
velocity_field, velocity_field, common_data, sit.second));
fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsDiagNonLin(
velocity_field, velocity_field, common_data, sit.second));
fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsOffDiag(
velocity_field, pressure_field, common_data, sit.second));
fe_rhs_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
velocity_field, common_data->velPtr));
fe_rhs_ptr->getOpPtrVector().push_back(
common_data->gradVelPtr));
fe_rhs_ptr->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
pressure_field, common_data->pressPtr));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpAssembleRhsVelocityLin(velocity_field, common_data, sit.second));
fe_rhs_ptr->getOpPtrVector().push_back(new OpAssembleRhsVelocityNonLin(
velocity_field, common_data, sit.second));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpAssembleRhsPressure(pressure_field, common_data, sit.second));
}
};
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_rhs_ptr,
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs_ptr,
const std::string velocity_field, const std::string pressure_field,
boost::shared_ptr<CommonData> common_data, const EntityType type) {
for (auto &sit : common_data->setOfBlocksData) {
if (type == MBPRISM) {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
fe_lhs_ptr->getOpPtrVector().push_back(
new OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms());
fe_lhs_ptr->getOpPtrVector().push_back(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
fe_lhs_ptr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms());
fe_rhs_ptr->getOpPtrVector().push_back(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
}
fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsDiagLin(
velocity_field, velocity_field, common_data, sit.second));
fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsOffDiag(
velocity_field, pressure_field, common_data, sit.second));
fe_rhs_ptr->getOpPtrVector().push_back(
common_data->gradVelPtr));
fe_rhs_ptr->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
pressure_field, common_data->pressPtr));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpAssembleRhsVelocityLin(velocity_field, common_data, sit.second));
fe_rhs_ptr->getOpPtrVector().push_back(
new OpAssembleRhsPressure(pressure_field, common_data, sit.second));
}
};
boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_flux_ptr,
const std::string velocity_field, boost::shared_ptr<CommonData> common_data,
const EntityType type) {
for (auto &sit : common_data->setOfBlocksData) {
if (type == MBPRISM) {
boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
fe_flux_ptr->getOpPtrVector().push_back(
new OpMultiplyDeterminantOfJacobianAndWeightsForFatPrisms());
fe_flux_ptr->getOpPtrVector().push_back(
new OpCalculateInvJacForFatPrism(inv_jac_ptr));
fe_flux_ptr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
}
fe_flux_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
velocity_field, common_data->velPtr));
fe_flux_ptr->getOpPtrVector().push_back(
new OpCalcVolumeFlux(velocity_field, common_data, sit.second));
}
};
boost::shared_ptr<FaceElementForcesAndSourcesCore> dragFe,
boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideDragFe,
std::string side_fe_name, const std::string velocity_field,
const std::string pressure_field,
boost::shared_ptr<CommonData> common_data) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
for (auto &sit : common_data->setOfFacesData) {
sideDragFe->getOpPtrVector().push_back(
common_data->gradVelPtr));
dragFe->getOpPtrVector().push_back(new OpCalculateHOJacForFace(jac_ptr));
dragFe->getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
dragFe->getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
dragFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
pressure_field, common_data->pressPtr));
dragFe->getOpPtrVector().push_back(
velocity_field, sideDragFe, side_fe_name, common_data, sit.second));
dragFe->getOpPtrVector().push_back(new NavierStokesElement::OpCalcDragForce(
velocity_field, common_data, sit.second));
}
};
boost::shared_ptr<PostProcFaceOnRefinedMesh> postProcDragPtr,
boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideDragFe,
std::string side_fe_name, const std::string velocity_field,
const std::string pressure_field,
boost::shared_ptr<CommonData> common_data) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
for (auto &sit : common_data->setOfFacesData) {
sideDragFe->getOpPtrVector().push_back(
common_data->gradVelPtr));
postProcDragPtr->getOpPtrVector().push_back(
new OpCalculateHOJacForFace(jac_ptr));
postProcDragPtr->getOpPtrVector().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
postProcDragPtr->getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(inv_jac_ptr));
postProcDragPtr->getOpPtrVector().push_back(
new OpCalculateVectorFieldValues<3>(velocity_field,
common_data->velPtr));
postProcDragPtr->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues(pressure_field,
common_data->pressPtr));
postProcDragPtr->getOpPtrVector().push_back(
velocity_field, sideDragFe, side_fe_name, common_data, sit.second));
postProcDragPtr->getOpPtrVector().push_back(
velocity_field, postProcDragPtr->postProcMesh,
postProcDragPtr->mapGaussPts, common_data, sit.second));
}
};
int row_side, int col_side, EntityType row_type, EntityType col_type,
EntData &row_data, EntData &col_data) {
row_nb_dofs = row_data.getIndices().size();
if (!row_nb_dofs)
col_nb_dofs = col_data.getIndices().size();
if (!col_nb_dofs)
if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
blockData.eNts.end()) {
}
row_nb_gauss_pts = row_data.getN().size1();
if (!row_nb_gauss_pts) {
}
isOnDiagonal = (row_type == col_type) && (row_side == col_side);
// Set size can clear local tangent matrix
locMat.resize(row_nb_dofs, col_nb_dofs, false);
locMat.clear();
CHKERR iNtegrate(row_data, col_data);
CHKERR aSsemble(row_data, col_data);
}
EntData &col_data) {
const int *row_ind = &*row_data.getIndices().begin();
const int *col_ind = &*col_data.getIndices().begin();
Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
: getFEMethod()->snes_B;
CHKERR MatSetValues(B, row_nb_dofs, row_ind, col_nb_dofs, col_ind,
&*locMat.data().begin(), ADD_VALUES);
if (!diagonalBlock || (sYmm && !isOnDiagonal)) {
locMat = trans(locMat);
CHKERR MatSetValues(B, col_nb_dofs, col_ind, row_nb_dofs, row_ind,
&*locMat.data().begin(), ADD_VALUES);
}
}
EntData &col_data) {
const int row_nb_base_functions = row_data.getN().size2();
auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
// INTEGRATION
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
int row_bb = 0;
for (; row_bb != row_nb_dofs / 3; row_bb++) {
t1(i) = w * row_diff_base_functions(i);
auto base_functions = col_data.getFTensor0N(gg, 0);
for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
FTensor::Tensor1<double *, 3> k(&locMat(3 * row_bb + 0, col_bb),
&locMat(3 * row_bb + 1, col_bb),
&locMat(3 * row_bb + 2, col_bb));
k(i) -= t1(i) * base_functions;
++base_functions;
}
++row_diff_base_functions;
}
for (; row_bb != row_nb_base_functions; row_bb++) {
++row_diff_base_functions;
}
}
}
EntData &col_data) {
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
&m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
&m(r + 2, c + 2));
};
const int row_nb_base_functions = row_data.getN().size2();
auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
// integrate local matrix for entity block
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
double const alpha = w * blockData.viscousCoef;
int row_bb = 0;
for (; row_bb != row_nb_dofs / 3; row_bb++) {
auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
const int final_bb = isOnDiagonal ? row_bb + 1 : col_nb_dofs / 3;
int col_bb = 0;
for (; col_bb != final_bb; col_bb++) {
auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
for (int i : {0, 1, 2}) {
for (int j : {0, 1, 2}) {
t_assemble(i, i) +=
alpha * row_diff_base_functions(j) * col_diff_base_functions(j);
t_assemble(i, j) +=
alpha * row_diff_base_functions(j) * col_diff_base_functions(i);
}
}
// Next base function for column
++col_diff_base_functions;
}
// Next base function for row
++row_diff_base_functions;
}
for (; row_bb != row_nb_base_functions; row_bb++) {
++row_diff_base_functions;
}
}
if (isOnDiagonal) {
for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
int col_bb = 0;
for (; col_bb != row_bb + 1; col_bb++) {
auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
auto t_off_side = get_tensor2(locMat, 3 * col_bb, 3 * row_bb);
t_off_side(i, j) = t_assemble(j, i);
}
}
}
}
EntData &col_data) {
auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
&m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
&m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
&m(r + 2, c + 2));
};
const int row_nb_base_functions = row_data.getN().size2();
auto row_base_functions = row_data.getFTensor0N();
auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
// integrate local matrix for entity block
for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
// Get volume and integration weight
double w = getVolume() * getGaussPts()(3, gg);
double const beta = w * blockData.inertiaCoef;
int row_bb = 0;
for (; row_bb != row_nb_dofs / 3; row_bb++) {
auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
auto col_base_functions = col_data.getFTensor0N(gg, 0);
const int final_bb = col_nb_dofs / 3;
int col_bb = 0;
for (; col_bb != final_bb; col_bb++) {
auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
for (int i : {0, 1, 2}) {
for (int j : {0, 1, 2}) {
t_assemble(i, j) +=
beta * col_base_functions * t_u_grad(i, j) * row_base_functions;
t_assemble(i, i) +=
beta * t_u(j) * col_diff_base_functions(j) * row_base_functions;
}
}
// Next base function for column
++col_diff_base_functions;
++col_base_functions;
}
// Next base function for row
++row_base_functions;
}
for (; row_bb != row_nb_base_functions; row_bb++) {
++row_base_functions;
}
++t_u;
++t_u_grad;
}
}
EntityType row_type,
EntData &row_data) {
// get number of dofs on row
nbRows = row_data.getIndices().size();
if (!nbRows)
if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
blockData.eNts.end()) {
}
// get number of integration points
nbIntegrationPts = getGaussPts().size2();
// integrate local vector
CHKERR iNtegrate(row_data);
// assemble local vector
CHKERR aSsemble(row_data);
}
// get global indices of local vector
const int *indices = &*data.getIndices().data().begin();
// get values from local vector
const double *vals = &*locVec.data().begin();
Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
: getFEMethod()->snes_f;
// assemble vector
CHKERR VecSetValues(f, nbRows, indices, vals, ADD_VALUES);
}
auto get_tensor1 = [](VectorDouble &v, const int r) {
&v(r + 0), &v(r + 1), &v(r + 2));
};
// set size of local vector
locVec.resize(nbRows, false);
// clear local entity vector
locVec.clear();
int nb_base_functions = data.getN().size2();
// get base function gradient on rows
auto t_v_grad = data.getFTensor1DiffN<3>();
auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
auto t_p = getFTensor0FromVec(*commonData->pressPtr);
// loop over all integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
double w = getVolume() * getGaussPts()(3, gg);
// evaluate constant term
const double alpha = w * blockData.viscousCoef;
auto t_a = get_tensor1(locVec, 0);
int rr = 0;
// loop over base functions
for (; rr != nbRows / 3; rr++) {
// add to local vector source term
t_a(i) += alpha * t_u_grad(i, j) * t_v_grad(j);
t_a(j) += alpha * t_u_grad(i, j) * t_v_grad(i);
t_a(i) -= w * t_p * t_v_grad(i);
++t_a; // move to next element of local vector
++t_v_grad; // move to next gradient of base function
}
for (; rr != nb_base_functions; rr++) {
++t_v_grad;
}
++t_u_grad;
++t_p;
}
}
auto get_tensor1 = [](VectorDouble &v, const int r) {
&v(r + 0), &v(r + 1), &v(r + 2));
};
// set size of local vector
locVec.resize(nbRows, false);
// clear local entity vector
locVec.clear();
// get base functions on entity
auto t_v = data.getFTensor0N();
int nb_base_functions = data.getN().size2();
auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
// loop over all integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
double w = getVolume() * getGaussPts()(3, gg);
// evaluate constant term
const double beta = w * blockData.inertiaCoef;
auto t_a = get_tensor1(locVec, 0);
int rr = 0;
// loop over base functions
for (; rr != nbRows / 3; rr++) {
// add to local vector source term
t_a(i) += beta * t_v * t_u_grad(i, j) * t_u(j);
++t_a; // move to next element of local vector
++t_v; // move to next base function
}
for (; rr != nb_base_functions; rr++) {
++t_v;
}
++t_u;
++t_u_grad;
}
}
// commonData->getBlockData(blockData);
// set size to local vector
locVec.resize(nbRows, false);
// clear local vector
locVec.clear();
// get base function
auto t_q = data.getFTensor0N();
int nb_base_functions = data.getN().size2();
// get solution at integration point
auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
// FTensor::Index<'i', 3> i;
// make loop over integration points
for (int gg = 0; gg != nbIntegrationPts; gg++) {
// evaluate function on boundary and scale it by area and integration
// weight
double w = getVolume() * getGaussPts()(3, gg);
// get element of vector
int rr = 0;
for (; rr != nbRows; rr++) {
for (int ii : {0, 1, 2}) {
t_a -= w * t_q * t_u_grad(ii, ii);
}
++t_a;
++t_q;
}
for (; rr != nb_base_functions; rr++) {
++t_q;
}
++t_u_grad;
}
}
EntityType type,
EntData &data) {
if (type != MBVERTEX)
PetscFunctionReturn(0);
if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
blockData.eNts.end()) {
}
const int nb_gauss_pts = getGaussPts().size2();
auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
auto t_w = getFTensor0IntegrationWeight(); ///< Integration weight
t_flux(i) = 0.0; // Zero entries
// loop over all integration points
for (int gg = 0; gg != nb_gauss_pts; ++gg) {
double vol = getVolume();
t_flux(i) += t_w * vol * t_u(i);
++t_w;
++t_u;
}
// Set array of indices
constexpr std::array<int, 3> indices = {0, 1, 2};
// Assemble volumetric flux
CHKERR VecSetValues(commonData->volumeFluxVec, 3, indices.data(), &t_flux(0),
ADD_VALUES);
}
EntityType type,
EntData &data) {
if (type != MBVERTEX)
PetscFunctionReturn(0);
if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
blockData.eNts.end()) {
}
const int nb_gauss_pts = getGaussPts().size2();
auto pressure_drag_at_gauss_pts =
getFTensor1FromMat<3>(*commonData->pressureDragTract);
auto shear_drag_at_gauss_pts =
getFTensor1FromMat<3>(*commonData->shearDragTract);
auto total_drag_at_gauss_pts =
getFTensor1FromMat<3>(*commonData->totalDragTract);
FTensor::Tensor1<double, 3> t_pressure_drag_force;
FTensor::Tensor1<double, 3> t_shear_drag_force;
FTensor::Tensor1<double, 3> t_total_drag_force;
t_pressure_drag_force(i) = 0.0; // Zero entries
t_shear_drag_force(i) = 0.0; // Zero entries
t_total_drag_force(i) = 0.0; // Zero entries
auto t_w = getFTensor0IntegrationWeight();
auto t_normal = getFTensor1NormalsAtGaussPts();
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double nrm2 = std::sqrt(t_normal(i) * t_normal(i));
double w = t_w * nrm2 * 0.5;
t_pressure_drag_force(i) += w * pressure_drag_at_gauss_pts(i);
t_shear_drag_force(i) += w * shear_drag_at_gauss_pts(i);
t_total_drag_force(i) += w * total_drag_at_gauss_pts(i);
++t_w;
++t_normal;
++pressure_drag_at_gauss_pts;
++shear_drag_at_gauss_pts;
++total_drag_at_gauss_pts;
}
// Set array of indices
constexpr std::array<int, 3> indices = {0, 1, 2};
CHKERR VecSetValues(commonData->pressureDragForceVec, 3, indices.data(),
&t_pressure_drag_force(0), ADD_VALUES);
CHKERR VecSetValues(commonData->shearDragForceVec, 3, indices.data(),
&t_shear_drag_force(0), ADD_VALUES);
CHKERR VecSetValues(commonData->totalDragForceVec, 3, indices.data(),
&t_total_drag_force(0), ADD_VALUES);
}
EntityType type,
EntData &data) {
if (type != MBVERTEX)
PetscFunctionReturn(0);
if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
blockData.eNts.end()) {
}
CHKERR loopSideVolumes(sideFeName, *sideFe);
const int nb_gauss_pts = getGaussPts().size2();
auto t_p = getFTensor0FromVec(*commonData->pressPtr);
auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
commonData->pressureDragTract->resize(3, nb_gauss_pts, false);
commonData->pressureDragTract->clear();
commonData->shearDragTract->resize(3, nb_gauss_pts, false);
commonData->shearDragTract->clear();
commonData->totalDragTract->resize(3, nb_gauss_pts, false);
commonData->totalDragTract->clear();
auto pressure_drag_at_gauss_pts =
getFTensor1FromMat<3>(*commonData->pressureDragTract);
auto shear_drag_at_gauss_pts =
getFTensor1FromMat<3>(*commonData->shearDragTract);
auto total_drag_at_gauss_pts =
getFTensor1FromMat<3>(*commonData->totalDragTract);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
double nrm2 = std::sqrt(t_normal(i) * t_normal(i));
t_unit_normal(i) = t_normal(i) / nrm2;
double mu = blockData.fluidViscosity;
pressure_drag_at_gauss_pts(i) = t_p * t_unit_normal(i);
shear_drag_at_gauss_pts(i) =
-mu * (t_u_grad(i, j) + t_u_grad(j, i)) * t_unit_normal(j);
total_drag_at_gauss_pts(i) =
pressure_drag_at_gauss_pts(i) + shear_drag_at_gauss_pts(i);
++pressure_drag_at_gauss_pts;
++shear_drag_at_gauss_pts;
++total_drag_at_gauss_pts;
++t_p;
++t_u_grad;
++t_normal;
}
}
EntityType type,
EntData &data) {
if (type != MBVERTEX)
PetscFunctionReturn(0);
double def_VAL[9];
bzero(def_VAL, 9 * sizeof(double));
if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
blockData.eNts.end()) {
}
Tag th_velocity;
Tag th_pressure;
Tag th_velocity_grad;
Tag th_shear_drag;
Tag th_pressure_drag;
Tag th_total_drag;
CHKERR postProcMesh.tag_get_handle("VELOCITY", 3, MB_TYPE_DOUBLE, th_velocity,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle("PRESSURE", 1, MB_TYPE_DOUBLE, th_pressure,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle("VELOCITY_GRAD", 9, MB_TYPE_DOUBLE,
th_velocity_grad,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle("PRESSURE_DRAG", 3, MB_TYPE_DOUBLE,
th_pressure_drag,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle("SHEAR_DRAG", 3, MB_TYPE_DOUBLE,
th_shear_drag,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle("TOTAL_DRAG", 3, MB_TYPE_DOUBLE,
th_total_drag,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
auto t_p = getFTensor0FromVec(*commonData->pressPtr);
const int nb_gauss_pts = commonData->pressureDragTract->size2();
MatrixDouble velGradMat;
velGradMat.resize(3, 3);
VectorDouble velVec;
velVec.resize(3);
VectorDouble pressDragVec;
pressDragVec.resize(3);
VectorDouble viscDragVec;
viscDragVec.resize(3);
VectorDouble totDragVec;
totDragVec.resize(3);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
for (int i : {0, 1, 2}) {
for (int j : {0, 1, 2}) {
velGradMat(i, j) = (*commonData->gradVelPtr)(i * 3 + j, gg);
}
velVec(i) = (*commonData->velPtr)(i, gg);
pressDragVec(i) = (*commonData->pressureDragTract)(i, gg);
viscDragVec(i) = (*commonData->shearDragTract)(i, gg);
totDragVec(i) = (*commonData->totalDragTract)(i, gg);
}
CHKERR postProcMesh.tag_set_data(th_velocity, &mapGaussPts[gg], 1,
&velVec(0));
CHKERR postProcMesh.tag_set_data(th_pressure, &mapGaussPts[gg], 1, &t_p);
CHKERR postProcMesh.tag_set_data(th_velocity_grad, &mapGaussPts[gg], 1,
&velGradMat(0, 0));
CHKERR postProcMesh.tag_set_data(th_pressure_drag, &mapGaussPts[gg], 1,
&pressDragVec(0));
CHKERR postProcMesh.tag_set_data(th_shear_drag, &mapGaussPts[gg], 1,
&viscDragVec(0));
CHKERR postProcMesh.tag_set_data(th_total_drag, &mapGaussPts[gg], 1,
&totDragVec(0));
++t_p;
}
}
EntityType type,
EntData &data) {
if (type != MBVERTEX)
PetscFunctionReturn(0);
double def_VAL[9];
bzero(def_VAL, 9 * sizeof(double));
if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
blockData.eNts.end()) {
}
// commonData->getBlockData(blockData);
Tag th_vorticity;
Tag th_q;
Tag th_l2;
CHKERR postProcMesh.tag_get_handle("VORTICITY", 3, MB_TYPE_DOUBLE,
th_vorticity, MB_TAG_CREAT | MB_TAG_SPARSE,
def_VAL);
CHKERR postProcMesh.tag_get_handle("Q", 1, MB_TYPE_DOUBLE, th_q,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
CHKERR postProcMesh.tag_get_handle("L2", 1, MB_TYPE_DOUBLE, th_l2,
MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
// auto p = getFTensor0FromVec(*commonData->pressPtr);
const int nb_gauss_pts = commonData->gradVelPtr->size2();
// const int nb_gauss_pts2 = commonData->pressPtr->size();
// const double lambda = commonData->lAmbda;
// const double mu = commonData->mU;
// FTensor::Index<'i', 3> i;
// FTensor::Index<'j', 3> j;
// FTensor::Index<'j', 3> k;
// FTensor::Tensor2<double,3,3> t_s;
double q;
double l2;
// FTensor::Tensor2<double, 3, 3> stress;
MatrixDouble Omega;
S.resize(3, 3);
Omega.resize(3, 3);
M.resize(3, 3);
for (int gg = 0; gg != nb_gauss_pts; gg++) {
vorticity(0) = t_u_grad(2, 1) - t_u_grad(1, 2);
vorticity(1) = t_u_grad(0, 2) - t_u_grad(2, 0);
vorticity(2) = t_u_grad(1, 0) - t_u_grad(0, 1);
q = 0;
for (int i = 0; i != 3; i++) {
for (int j = 0; j != 3; j++) {
q += -0.5 * t_u_grad(i, j) * t_u_grad(j, i);
}
}
for (int i = 0; i != 3; i++) {
for (int j = 0; j != 3; j++) {
S(i, j) = 0.5 * (t_u_grad(i, j) + t_u_grad(j, i));
Omega(i, j) = 0.5 * (t_u_grad(i, j) - t_u_grad(j, i));
M(i, j) = 0.0;
}
}
for (int i = 0; i != 3; i++) {
for (int j = 0; j != 3; j++) {
for (int k = 0; k != 3; k++) {
M(i, j) += S(i, k) * S(k, j) + Omega(i, k) * Omega(k, j);
}
}
}
MatrixDouble eigen_vectors = M;
VectorDouble eigen_values(3);
// LAPACK - eigenvalues and vectors. Applied twice for initial creates
// memory space
int n = 3, lda = 3, info, lwork = -1;
double wkopt;
info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
&(eigen_values.data()[0]), &wkopt, lwork);
if (info != 0)
SETERRQ1(PETSC_COMM_SELF, 1,
"is something wrong with lapack_dsyev info = %d", info);
lwork = (int)wkopt;
double work[lwork];
info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
&(eigen_values.data()[0]), work, lwork);
if (info != 0)
SETERRQ1(PETSC_COMM_SELF, 1,
"is something wrong with lapack_dsyev info = %d", info);
map<double, int> eigen_sort;
eigen_sort[eigen_values[0]] = 0;
eigen_sort[eigen_values[1]] = 1;
eigen_sort[eigen_values[2]] = 2;
// prin_stress_vect.clear();
VectorDouble prin_vals_vect(3);
prin_vals_vect.clear();
int ii = 0;
for (map<double, int>::reverse_iterator mit = eigen_sort.rbegin();
mit != eigen_sort.rend(); mit++) {
prin_vals_vect[ii] = eigen_values[mit->second];
// for (int dd = 0; dd != 3; dd++) {
// prin_stress_vect(ii, dd) = eigen_vectors.data()[3 * mit->second +
// dd];
// }
ii++;
}
l2 = prin_vals_vect[1];
// cout << prin_vals_vect << endl;
// cout << "-0.5 sum: " << -0.5 * (prin_vals_vect[0] + prin_vals_vect[1]
// + prin_vals_vect[2]) << endl; cout << "q: " << q << endl;
// t_s(i,j) = 0.5*(t)
// vorticity(0) = t_u_grad(1, 2) - t_u_grad(2, 1);
// vorticity(1) = t_u_grad(2, 0) - t_u_grad(0, 2);
// vorticity(2) = t_u_grad(0, 1) - t_u_grad(1, 0);
CHKERR postProcMesh.tag_set_data(th_vorticity, &mapGaussPts[gg], 1,
&vorticity(0));
CHKERR postProcMesh.tag_set_data(th_q, &mapGaussPts[gg], 1, &q);
CHKERR postProcMesh.tag_set_data(th_l2, &mapGaussPts[gg], 1, &l2);
++t_u_grad;
}
}
VectorDouble3 stokes_flow_velocity(double x, double y, double z) {
double r = std::sqrt(x * x + y * y + z * z);
double theta = acos(x / r);
double phi = atan2(y, z);
double ur = cos(theta) * (1.0 + 0.5 / (r * r * r) - 1.5 / r);
double ut = -sin(theta) * (1.0 - 0.25 / (r * r * r) - 0.75 / r);
VectorDouble3 res(3);
res[0] = ur * cos(theta) - ut * sin(theta);
res[1] = ur * sin(theta) * sin(phi) + ut * cos(theta) * sin(phi);
res[2] = ur * sin(theta) * cos(phi) + ut * cos(theta) * cos(phi);
return res;
}
static Index< 'M', 3 > M
static Index< 'n', 3 > n
VectorDouble3 stokes_flow_velocity(double x, double y, double z)
constexpr double alpha
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:273
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:103
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
double w(const double g, const double t)
const double r
rate factor
static double phi
constexpr double mu
FTensor::Index< 'm', 3 > m
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode aSsemble(EntData &data)
MoFEMErrorCode iNtegrate(EntData &data)
Integrate local constrains vector.
MoFEMErrorCode iNtegrate(EntData &data)
Integrate local entity vector.
Calculate drag force on the fluid-solid interface.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate drag traction on the fluid-solid interface.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Post processing output of drag traction on the fluid-solid interface.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
static MoFEMErrorCode setCalcDragOperators(boost::shared_ptr< FaceElementForcesAndSourcesCore > dragFe, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for calculating drag force on the solid surface.
static MoFEMErrorCode setPostProcDragOperators(boost::shared_ptr< PostProcFaceOnRefinedMesh > postProcDragPtr, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for post processing output of drag traction.
static MoFEMErrorCode setNavierStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Navier-Stokes equations.
static MoFEMErrorCode setStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Stokes equations.
static MoFEMErrorCode setCalcVolumeFluxOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe_flux_ptr, const std::string velocity_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for calculation of volume flux.