Commit 60438517 authored by Tomáš Jakubec's avatar Tomáš Jakubec
Browse files

mesh data export fix

parent 9c1dd0f2
......@@ -221,7 +221,7 @@ void testHeatConduction1() {
constexpr unsigned int ProblemDim = 3;
MultiphaseFlow mpf;
mpf.setupMeshData("cube_3500_cells.vtk");
mpf.setupMeshData("cube_64k_cells.vtk");
FlowData<ProblemDim> initGlobals;
initGlobals.R_spec = 287;
......@@ -246,7 +246,7 @@ void testHeatConduction1() {
mpf.inFlow_eps_g = 1;
mpf.inFlow_eps_s = 0;
mpf.inFlow_u_g = {0.0,25};
mpf.inFlow_u_g = {0.0,5};
mpf.inFlow_u_s = {0, 0};
......@@ -262,20 +262,21 @@ void testHeatConduction1() {
MeshDataContainer<FlowData<ProblemDim>, ProblemDim> compData(mpf.mesh, ini);
for (auto& cell : mpf.mesh.getCells()){
if(cell.getCenter()[1] > 1.0 && cell.getCenter()[1] < 7.0 &&
cell.getCenter()[0] > 0.0 && cell.getCenter()[0] < 10.0){
if(cell.getCenter()[2] > -0.5 && cell.getCenter()[2] < 0.5 &&
cell.getCenter()[1] > -0.5 && cell.getCenter()[1] < 0.5 &&
cell.getCenter()[0] > -0.5 && cell.getCenter()[0] < 0.5){
compData.at(cell).eps_s = 0.2;
}
}
DBGVAR_JSON(mpf.mesh.computeElementMeasures().getDataByDim<3>());
mpf.exportData(0.0, compData);
return;
double exportStep = 1e-1;
for (double t = 0; t < /*300 * */exportStep; t += exportStep){
for (double t = 0; t < 10 * exportStep; t += exportStep){
RKMSolver(mpf, compData, 1e-3, t, t + exportStep, 1e-1);
RKMSolver(mpf, compData, 1e-3, t, t + exportStep, 1e-2);
//EulerSolver(mpf, compData, 1e-5, t, t + exportStep);
mpf.exportData((t + exportStep)*1e-2, compData);
......
......@@ -256,7 +256,9 @@ void MultiphaseFlow::ComputeSource(const MeshType::Cell& cell,
Vector<ProblemDimension, double> drag = Beta_s(cellData) * (cellData.getVelocitySolid() - cellData.getVelocityGas());
Vector<ProblemDimension,double> g_acceleration = {0, -9.81};
Vector<ProblemDimension,double> g_acceleration = {};
g_acceleration[g_acceleration.size() - 1] = -9.81;
resData.p_g += (cellData.rho_g * g_acceleration + drag);
......@@ -280,10 +282,11 @@ void MultiphaseFlow::ComputeSource(const MeshType::Cell& cell,
double MultiphaseFlow::FlowModulation(const Vertex<MeshType::meshDimension(), double>& x)
{
double inFlowModulation = x[0];
double inFlowModulation = sqrt(pow(x[0],2) + pow(x[1],2));
//inFlowModulation = (- inFlowModulation * inFlowModulation + inFlowModulation - 0.0099) * 4.164931279;
inFlowModulation = -(inFlowModulation -1.9) * (inFlowModulation - 4.7) * 0.5102;//(-inFlowModulation * inFlowModulation + 6.6 * inFlowModulation - 8.93) * 0.5102;
//inFlowModulation = -(inFlowModulation -1.9) * (inFlowModulation - 4.7) * 0.5102;//(-inFlowModulation * inFlowModulation + 6.6 * inFlowModulation - 8.93) * 0.5102;
inFlowModulation = (inFlowModulation - 0.3) * (inFlowModulation - 0.3) * (1/0.09);
return inFlowModulation;
}
......@@ -381,19 +384,24 @@ void MultiphaseFlow::setupMeshData(const std::string& fileName){
Type MultiphaseFlow::TypeOfCell(const MeshType::Cell &cell) {
if (cell.getIndex() >= BOUNDARY_INDEX(size_t)){
/*
if (cell.getCenter()[1] <= 1e-5) {
if (
sqrt(pow(cell.getCenter()[0],2) + pow(cell.getCenter()[1],2)) < 0.3 &&
cell.getCenter()[2] < 0
) {
return Type::INFLOW;
}
if (cell.getCenter()[1] >= 34.349
if (
sqrt(pow(cell.getCenter()[0],2) + pow(cell.getCenter()[1],2)) < 0.3 &&
cell.getCenter()[2] > 0
//cell.GetCenter()[0] > 2 && (cell.GetCenter()[1] < 32.0 && cell.GetCenter()[1] > 30)
//cell.GetCenter()[0] > 8 && cell.GetCenter()[1] >= 33.39
//cell.getCenter()[1] >= 1.999
) {
return Type::OUTFLOW;
}
*/
return Type::WALL;
} else {
......
......@@ -649,12 +649,12 @@ inline void MultiphaseFlow::ComputeFluxGas_inner(const FlowData<ProblemDimension
}
inline void MultiphaseFlow::ComputeFluxGas_inflow(const FlowData<ProblemDimension>& innerCellData, const Cell& innerCell, FaceData<ProblemDimension>& edgeData, const Face&)
inline void MultiphaseFlow::ComputeFluxGas_inflow(const FlowData<ProblemDimension>& innerCellData, const Cell& innerCell, FaceData<ProblemDimension>& edgeData, const Face& face)
{
Vector<ProblemDimension,double> modulatedU = inFlow_u_g;
modulatedU[1] *= FlowModulation(innerCell.getCenter());
modulatedU[1] *= FlowModulation(face.getCenter());
double productOf_u_And_n = (modulatedU * edgeData.n);
......@@ -865,13 +865,13 @@ inline void MultiphaseFlow::ComputeFluxSolid_inner(const FlowData<ProblemDimensi
}
inline void MultiphaseFlow::ComputeFluxSolid_inflow(const FlowData<ProblemDimension>& innerCellData, const Cell& innerCell, FaceData<ProblemDimension>& edgeData, const Face&)
inline void MultiphaseFlow::ComputeFluxSolid_inflow(const FlowData<ProblemDimension>& innerCellData, const Cell& innerCell, FaceData<ProblemDimension>& edgeData, const Face& face)
{
(void)innerCellData;
Vector<ProblemDimension,double> modulatedU = inFlow_u_s;
modulatedU[1] *= FlowModulation(innerCell.getCenter());
modulatedU[1] *= FlowModulation(face.getCenter());
double productOf_u_And_n = (modulatedU * edgeData.n);
......
......@@ -40,7 +40,9 @@ class VTKMeshDataWriter {
IndexType realIndex = 0;
IndexType localIndex = 0;
for(const std::pair<IndexType, IndexType>& key : writer.backwardCellIndexMapping) {
auto keyIt = writer.backwardCellIndexMapping.cbegin();
while(keyIt != writer.backwardCellIndexMapping.cend()) {
const std::pair<IndexType, IndexType>& key = *keyIt;
while (localIndex < key.first) {
for (unsigned int j = 0; j < DefaultIOTraits<T>::getTraits().template getValue<Index>(data.at(0)).size(); j++) {
ost << DefaultIOTraits<T>::getTraits().template getValue<Index>(data.at(realIndex))[j] << ' ';
......@@ -53,7 +55,13 @@ class VTKMeshDataWriter {
for (unsigned int j = 0; j < DefaultIOTraits<T>::getTraits().template getValue<Index>(data.at(0)).size(); j++) {
ost << DefaultIOTraits<T>::getTraits().template getValue<Index>(data.at(realIndex))[j] << ' ';
}
++keyIt;
if (keyIt == writer.backwardCellIndexMapping.cend() || realIndex != keyIt->second){
realIndex++;
}
}
// write the rest of not tessellated cells
while (realIndex < data.size()) {
for (unsigned int j = 0; j < DefaultIOTraits<T>::getTraits().template getValue<Index>(data.at(0)).size(); j++) {
ost << DefaultIOTraits<T>::getTraits().template getValue<Index>(data.at(realIndex))[j] << ' ';
......@@ -75,7 +83,10 @@ class VTKMeshDataWriter {
IndexType realIndex = 0;
IndexType localIndex = 0;
for(const std::pair<IndexType, IndexType>& key : writer.backwardCellIndexMapping) {
auto keyIt = writer.backwardCellIndexMapping.cbegin();
while(keyIt != writer.backwardCellIndexMapping.cend()) {
const std::pair<IndexType, IndexType>& key = *keyIt;
while (localIndex < key.first) {
ost << DefaultIOTraits<T>::getTraits().template getValue<Index>(data.at(realIndex)) << ' ';
realIndex++;
......@@ -84,7 +95,12 @@ class VTKMeshDataWriter {
realIndex = key.second;
localIndex++;
ost << DefaultIOTraits<T>::getTraits().template getValue<Index>(data.at(realIndex)) << ' ';
++keyIt;
if (keyIt == writer.backwardCellIndexMapping.cend() || realIndex != keyIt->second){
realIndex++;
}
}
while (realIndex < data.size()) {
ost << DefaultIOTraits<T>::getTraits().template getValue<Index>(data.at(realIndex)) << ' ';
realIndex++;
......
......@@ -191,7 +191,6 @@ public:
}
};
template<>
class VTKMeshReader<3> : public MeshReader<3>{
using reader = MeshReader<3>;
......@@ -219,7 +218,7 @@ class VTKMeshReader<3> : public MeshReader<3>{
{4,0,5,8},
{5,1,6,9},
{6,2,7,10},
{7,3,5,11},
{7,3,4,11},
{8,9,10,11}
}
}
......@@ -375,9 +374,11 @@ public:
mesh.getFaces().at(faceIndex).setCellRightIndex(cellIndex);
}
// connecting of faces in cell structure
if (prevFaceIndex != INVALID_INDEX(IndexType)) {
mesh.getFaces().at(prevFaceIndex).setNextBElem(faceIndex, cellIndex);
}
if (fi == 0) {
mesh.getCells().at(cellIndex).setBoundaryElementIndex(faceIndex);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment