The computation is fixed and run in parallel. (132bcf77) · Commits · Tomáš Jakubec / GTMesh

Unstructured_mesh/Unstructured_mesh.pro

+2 −0

Original line number	Diff line number	Diff line
		@@ -2,6 +2,8 @@ TEMPLATE = app
		CONFIG += console c++17
		CONFIG -= app_bundle
		CONFIG -= qt
		QMAKE_CXXFLAGS += -fopenmp
		LIBS += -fopenmp

		SOURCES += \
		main.cpp \

Unstructured_mesh/main.cpp

+106 −197

Original line number	Diff line number	Diff line
		#include <iostream>
		//#define UNDEBUG
		#define CONSOLE_COLORED_OUTPUT
		//#define CONSOLE_COLORED_OUTPUT
		#include "../src/Debug/Debug.h"
		#include "multiphaseflow.h"

		@@ -12,129 +12,6 @@ using namespace std;



		// struct CellData {
		// double invVol;
		//
		// };
		//
		//
		// struct CompData {
		// double T;
		// //CompData(const CompData&) = default;
		// //CompData(CompData&&) = default;
		// CompData(double _T = 0){T = _T;}
		//
		// };
		// MAKE_NAMED_ATTRIBUTE_TRAIT(CompData, "Temperature", T);
		// MAKE_ATTRIBUTE_TRAIT_ARITHMETIC(CompData, T);
		//
		// struct FaceData {
		// double volOverDist;
		// };
		//
		// template <unsigned int ProblemDimension, typename Real>
		// class HeatCunduction {
		// public:
		// using MeshType = UnstructuredMesh<ProblemDimension,size_t, double, 12>;
		// using ResultType = CompData;
		//
		// public:
		// MeshType mesh;
		// MeshDataContainer<std::tuple<CellData, FaceData>, 3,2> meshData;
		// VTKMeshWriter<ProblemDimension, size_t, Real> writer;
		//
		// HeatCunduction(){
		// }
		//
		// void loadMesh(const std::string& fileName) {
		// FPMAMeshReader<MeshType::meshDimension()> reader;
		// ifstream file(fileName);
		// reader.loadFromStream(file, mesh);
		//
		// mesh.template initializeCenters<TESSELLATED>();
		//
		// mesh.setupBoundaryCells();
		//
		//
		//
		// mesh.setupBoundaryCellsCenters();
		//
		//
		// auto measures = mesh.template computeElementMeasures<TESSELLATED>();
		//
		// auto dists = ComputeCellsDistance(mesh);
		//
		//
		// meshData.allocateData(mesh);
		//
		// MeshDataContainer<CompData, 3> compData(mesh);
		//
		// for (auto& face : mesh.getFaces()) {
		// meshData.at(face).volOverDist = measures.at(face) / dists.at(face);
		// }
		//
		// for (auto& cell : mesh.getCells()) {
		// meshData.at(cell).invVol = 1.0 / measures.at(cell);
		// }
		// }
		//
		// void calculateRHS(
		// Real,//time is unused in this problem
		// MeshDataContainer<CompData, ProblemDimension>& compData,
		// MeshDataContainer<CompData, ProblemDimension>& outDeltas) {
		//
		// for (CompData& dT : outDeltas.template getDataByPos<0>()){
		// dT.T = 0;
		// }
		//
		// for (auto& face : mesh.getFaces()) {
		// CompData lT(1000);
		// bool lIn = false;
		// if(!((face.getCellLeftIndex() & BOUNDARY_INDEX(size_t)) == BOUNDARY_INDEX(size_t))){
		// lT = compData.template getDataByDim<3>().at(face.getCellLeftIndex());
		// lIn = true;
		// }
		//
		// CompData rT(1000);
		// bool rIn = false;
		// if(!((face.getCellRightIndex() & BOUNDARY_INDEX(size_t)) == BOUNDARY_INDEX(size_t))){
		// rT = compData.template getDataByDim<3>().at(face.getCellRightIndex());
		// rIn = true;
		// }
		//
		//
		// CompData dT = meshData.at(face).volOverDist * (lT - rT);
		// if (rIn)
		// outDeltas.template getDataByDim<3>().at(face.getCellRightIndex()) += dT;
		// if (lIn)
		// outDeltas.template getDataByDim<3>().at(face.getCellLeftIndex()) -= dT;
		//
		// }
		//
		// for (auto& cell : mesh.getCells()) {
		// outDeltas.at(cell) *= meshData.at(cell).invVol;
		// }
		// }
		//
		// template<typename dataType>
		// void exportData(double time,
		// MeshDataContainer<dataType, 3>& compData) {
		//
		// using std::string_literals;
		//
		// ofstream ofile("HeatCond"s + "_" + to_string(time) + ".vtk");
		// writer.writeHeader(ofile, "HC_test"s + to_string(time));
		// writer.writeToStream(ofile, mesh, MeshDataContainer<MeshNativeType<3>::ElementType, 3>(mesh, MeshNativeType<3>::POLYHEDRON));
		//
		// VTKMeshDataWriter<3> dataWriter;
		// dataWriter.writeToStream(ofile, compData, writer);
		//
		// ofile.close();
		// DBGMSG("Data eported");
		//
		// }
		// };

		struct watch{

		std::chrono::high_resolution_clock clock;
		@@ -210,43 +87,60 @@ void RKMSolver(
		MeshDataContainer<typename Problem::ResultType, MeshDimension> K5(compData);
		double time = startTime;
		bool run = true;

		DBGMSG("RKM_start", run);

		while (run) {
		#pragma omp single
		{
		wK1.start();
		if (time + tau > finalT) {
		tau = finalT - time;
		run = false;
		}
		}
		#pragma omp parallel
		{
		problem.calculateRHS(time, compData, K1);
		wK1.start();

		#pragma omp for
		for (size_t i = 0; i < Ktemp.template getDataByPos<0>().size(); i++){
		Ktemp.template getDataByPos<0>().at(i) = compData.template getDataByDim<MeshDimension>().at(i) + (tau * (1.0 / 3.0) * K1.template getDataByDim<MeshDimension>().at(i));
		}
		wK1.lap();


		problem.calculateRHS(time, Ktemp, K2);
		wK2.start();

		#pragma omp for
		for (size_t i = 0; i < Ktemp.template getDataByPos<0>().size(); i++){
		Ktemp.template getDataByPos<0>().at(i) = compData.template getDataByDim<MeshDimension>().at(i) + (tau * (1.0 / 6.0) * (K1.template getDataByDim<MeshDimension>().at(i) + K2.template getDataByDim<MeshDimension>().at(i)));
		}
		wK2.lap();


		problem.calculateRHS(time, Ktemp, K3);
		wK3.start();


		#pragma omp for
		for (size_t i = 0; i < Ktemp.template getDataByPos<0>().size(); i++){
		Ktemp.template getDataByPos<0>().at(i) = compData.template getDataByDim<MeshDimension>().at(i) + (tau * (0.125 * K1.template getDataByDim<MeshDimension>().at(i) + 0.375 * K3.template getDataByDim<MeshDimension>().at(i)));
		}
		wK3.lap();


		problem.calculateRHS(time, Ktemp, K4);
		wK4.start();

		#pragma omp for
		for (size_t i = 0; i < Ktemp.template getDataByPos<0>().size(); i++){
		Ktemp.template getDataByPos<0>().at(i) = compData.template getDataByDim<MeshDimension>().at(i) + (tau * ((0.5 * K1.template getDataByDim<MeshDimension>().at(i)) - (1.5 * K3.template getDataByDim<MeshDimension>().at(i)) + (2.0 * K4.template getDataByPos<0>().at(i))));
		}

		wK4.lap();
		problem.calculateRHS(time, Ktemp, K5);

		problem.calculateRHS(time, Ktemp, K5);
		}// end of parallel section

		double error = 0.0;
		wError.start();

		#pragma omp parallel
		#pragma omp for reduction(max : error)
		for (size_t i = 0; i < K4.template getDataByPos<0>().size(); i++){
		double tmpE = max(abs(0.2 * K1.template getDataByPos<0>().at(i) - 0.9 * K3.template getDataByPos<0>().at(i) +
		0.8 * K4.template getDataByPos<0>().at(i) - 0.1 * K5.template getDataByPos<0>().at(i)));
		@@ -255,33 +149,73 @@ wError.start();
		error = tmpE;
		}
		}
		wError.lap();



		error = tau (1.0 / 3.0);

		if (error < delta) {
		#pragma omp parallel
		#pragma omp for
		for (size_t i = 0; i < K4.template getDataByPos<0>().size(); i++){
		auto& _compData = compData.template getDataByPos<0>().at(i);
		_compData += tau * (1.0 / 6.0) * (((K1.template getDataByDim<MeshDimension>().at(i) + K5.template getDataByDim<MeshDimension>().at(i))) + (4.0 * K4.template getDataByPos<0>().at(i)));
		}

		time += tau;
		tau *= 1.005;
		//cout << "time: " << time << "\r";
		wK1.lap();

		} else {

		wK1.lap();
		tau = std::pow(0.2, delta/error) 0.8;
		run = true;

		}

		}

		DBGMSG("computation done");
		}


		template <typename Problem, typename = typename std::enable_if<HasDefaultArithmeticTraits<typename Problem::ResultType>::value>::type>
		void EulerSolver(
		Problem& problem,
		MeshDataContainer<typename Problem::ResultType, Problem::MeshType::meshDimension()>& compData,//x_ini
		double tau_ini,//tau_ini
		double startTime,
		double finalT
		)
		{
		constexpr unsigned int MeshDimension = Problem::MeshType::meshDimension();
		double tau = tau_ini;


		MeshDataContainer<typename Problem::ResultType, MeshDimension> K1(compData);
		double time = startTime;
		bool run = true;
		while (run) {
		if (time + tau >= finalT) {
		tau = finalT - time;
		run = false;
		}
		problem.calculateRHS(time, compData, K1);


		for (size_t i = 0; i < compData.template getDataByPos<0>().size(); i++){

		compData.template getDataByPos<0>().at(i) += tau * K1.template getDataByDim<MeshDimension>().at(i);
		}

		time += tau;
		}
		DBGMSG("computation done");
		}


		MAKE_ATTRIBUTE_TRAIT(CellData, invVolume);
		MAKE_ATTRIBUTE_TRAIT(EdgeData, n,LengthOverDist, Length,LeftCellKoef, RightCellKoef);
		MAKE_ATTRIBUTE_TRAIT(PointData, u_s, u_g, PointType, cellKoef);

		void testHeatConduction1() {

		@@ -293,97 +227,72 @@ void testHeatConduction1() {
		FlowData::T = 300;
		FlowData::rho_s = 1700;

		mpf.myu = 1e5;
		mpf.myu_s = 1.5;
		mpf.d_s = 0.06;
		mpf.outFlow.setPressure(1e5);
		//mpf.outFlow.rho_g = 1.3;
		//mpf.outFlow.eps_g = 1;
		mpf.outFlow.eps_s = 0;
		mpf.artifitialDisspation = 0.8;
		mpf.R_spec = 287;
		mpf.T = 300;
		mpf.artifitialDisspation = 1;
		mpf.phi_s = 1;
		mpf.myu = 1e-5;
		mpf.rho_s = 1700;
		mpf.myu_s = 0.5;//1.5;
		mpf.d_s = 0.002;//0.06;
		mpf.phi_s = 1;

		mpf.T = 300;




		mpf.outFlow.setPressure(1e5);
		mpf.outFlow.eps_s = 0;

		mpf.inFlow_eps_g = 1;
		mpf.inFlow_eps_s = 0;
		mpf.inFlow_u_g = {0.0,15};
		mpf.inFlow_u_g = {0.0,25};
		mpf.inFlow_u_s = {0, 0};


		FlowData ini;
		// ini.eps_g = 1;
		ini.eps_s = 0;
		ini.rho_g = 1.3;
		// ini.rho_g = 1.3;
		ini.setPressure(1e5);
		ini.p_g = {};
		//ini.setVelocityGas({0, 0});
		ini.p_s = {0, 0};


		MeshDataContainer<FlowData, 2> compData(mpf.mesh, ini);

		for (auto& cell : mpf.mesh.getCells()){
		if(cell.getCenter()[1] > 7 && cell.getCenter()[1] < 9){
		compData.at(cell).eps_s = 0.1;
		if(cell.getCenter()[1] > 1.0 && cell.getCenter()[1] < 7.0 &&
		cell.getCenter()[0] > 0.0 && cell.getCenter()[0] < 10.0){
		compData.at(cell).eps_s = 0.2;
		}
		}
		//MeshDataContainer<FlowData, 2> result(compData);

		//mpf.ActualizePointData(compData);

		//mpf.calculateRHS(0.0, compData, result);

		mpf.exportData(0.0, compData);
		for (double t = 0; t < 1; t += 0.05){
		RKMSolver(mpf, compData, 1e-4, t, t + 0.05, 1);
		mpf.exportData(t + 0.05, compData);
		}
		double exportStep = 1e-1;
		for (double t = 0; t < 300 * exportStep; t += exportStep){

		/*
		HeatCunduction<3,double> hcProblem;
		hcProblem.loadMesh("Poly_2_5_level2.fpma");

		double startTime = 0.0, finalTime = 0.1 , tau = 1e-4;
		RKMSolver(mpf, compData, 1e-3, t, t + exportStep, 1e-1);

		VTKMeshWriter<3, size_t, double> writer;

		MeshDataContainer<CompData, 3> compData(hcProblem.mesh, CompData(300));

		// hcProblem.exportData(startTime, compData);

		RKMSolver(hcProblem, compData,
		tau,
		startTime,
		finalTime,
		1);
		DBGVAR(wK1.getResult(),wK2.getResult(),wK3.getResult(),wK4.getResult(),wError.getResult());
		wK1.reset();wK2.reset();wK3.reset();wK4.reset();wError.reset();

		// hcProblem.exportData(finalTime, compData);



		startTime = finalTime;
		finalTime *= 2;
		RKMSolver(hcProblem, compData,
		tau,
		startTime,
		finalTime,
		1);
		//EulerSolver(mpf, compData, 1e-5, t, t + exportStep);
		mpf.exportData((t + exportStep)*1e-2, compData);

		DBGVAR(t + exportStep, wK1.getResult());
		}

		// hcProblem.exportData(finalTime, compData);

		DBGVAR(wK1.getResult(),wK2.getResult(),wK3.getResult(),wK4.getResult(),wError.getResult());
		*/


		}


		void atEnd(){
		DBGVAR(wK1.getResult());
		}

		int main()
		{
		atexit(atEnd);
		testHeatConduction1();
		//testHeatConduction();
		//meshAdmisibility("Poly_2_5_level1.fpma");

Unstructured_mesh/multiphaseflow.cpp

+41 −89

Original line number	Diff line number	Diff line
		@@ -14,53 +14,45 @@ void MultiphaseFlow::calculateRHS(double, MeshDataContainer<MultiphaseFlow::Resu

		//#pragma omp parallel
		{
		for (auto& cell : mesh.getCells()){

		#pragma omp for
		for (size_t cellIndex = 0; cellIndex < mesh.getCells().size(); cellIndex++){
		auto& cell = mesh.getCells()[cellIndex];
		FlowData& cellData = compData.at(cell);
		if (cellData.eps_s <= 0.0) {
		cellData.eps_s = 0.0;
		cellData.p_s = {};
		}
		}
		// first: compute vertexes values
		ActualizePointData(compData);
		// first: compute vertices values

		UpdateVertexData(compData);

		// second: compute fluxes over edges
		for (auto& face : mesh.getFaces())
		DBGTRY(ComputeFlux(face, compData);)
		//FluxComputationParallel(*this);
		//#pragma omp barrier
		//FVMethod::CellSourceComputationParallel(*this);
		for (auto& cell : mesh.getCells()){
		DBGTRY(ComputeSource(cell, compData, outDeltas);)
		}
		//#pragma omp barrier
		// third: compute values of new time step
		//FVMethod::TimeStepParallel(*this);
		}

		#pragma omp for
		for (size_t faceIndex = 0; faceIndex < mesh.getFaces().size(); faceIndex++){
		auto& face = mesh.getFaces()[faceIndex];

		ComputeFlux(face, compData);

		}

		// third: compute source terms

		/*
		void MultiphaseFlow::ComputeStep()
		{
		// in single compute step
		// we have to:
		#pragma omp for
		for (size_t cellIndex = 0; cellIndex < mesh.getCells().size(); cellIndex++){
		auto& cell = mesh.getCells()[cellIndex];

		//#pragma omp parallel
		{
		// first: compute vertexes values
		ActualizePointData();
		// second: compute fluxes over edges
		FluxComputationParallel(*this);
		//#pragma omp barrier
		FVMethod::CellSourceComputationParallel(*this);
		//#pragma omp barrier
		// third: compute values of new time step
		FVMethod::TimeStepParallel(*this);
		ComputeSource(cell, compData, outDeltas);
		}

		}
		*/

		}






		@@ -79,7 +71,6 @@ void MultiphaseFlow::ComputeFlux(const MeshType::Face &fcData, const MeshDataCon
		// Computation of fluxes of density and momentum of gaseous part
		ComputeFluxGas_inner(leftData, rightData, currEdgeData, fcData);


		// Compute flux of solid phase
		ComputeFluxSolid_inner(leftData, rightData, currEdgeData, fcData);

		@@ -156,7 +147,8 @@ void MultiphaseFlow::ComputeSource(const MeshType::Cell& ccData,
		mesh,
		// aplication of sum lambda to all cell faces
		[&](size_t cellIndex, size_t faceIndex){
		EdgeData& eData = meshData.getDataByDim<ProblemDimension - 1>().at(faceIndex);
		const EdgeData& eData = meshData.getDataByDim<ProblemDimension - 1>().at(faceIndex);

		if (cellIndex == mesh.getFaces().at(faceIndex).getCellLeftIndex()){
		cellData.fluxRho_g += eData.fluxRho_g;
		cellData.fluxRho_s += eData.fluxRho_s;
		@@ -171,13 +163,13 @@ void MultiphaseFlow::ComputeSource(const MeshType::Cell& ccData,
		}
		);


		cellData.fluxP_g *= meshData.at(ccData).invVolume;
		cellData.fluxRho_g *= meshData.at(ccData).invVolume;
		cellData.fluxP_s *= meshData.at(ccData).invVolume;
		cellData.fluxRho_s *= meshData.at(ccData).invVolume;



		Vector<ProblemDimension, double> drag = Beta_s(cellData) * (cellData.getVelocitySolid() - cellData.getVelocityGas());

		Vector<ProblemDimension,double> g_acceleration = {0, -9.81};
		@@ -212,53 +204,6 @@ void MultiphaseFlow::ComputeSource(const MeshType::Cell& ccData,



		// void MultiphaseFlow::ComputeTimeStep(FVMethod::CellComputeData &ccData)
		// {
		//
		// if (ccData.GetCurCell().GetCellFlag() == Type::WALL)
		// return;
		//
		//
		// FlowData& cellData = PrevState.GetDataAt(ccData.GetCurCellIndex());
		//
		//
		// // solid component
		// cellData.p_s += (cellData.fluxP_s * TimeStep);
		//
		// cellData.eps_s += cellData.fluxRho_s * (TimeStep / rho_s);
		//
		//
		//
		// /*
		// cellData.u_s = cellData.u_s + cellData.fluxP_s *
		// (TimeStep / reg(rho_s * cellData.eps_s)) ;
		// */
		// if (cellData.eps_s <= 0.0) {
		// cellData.eps_s = 0.0;
		// cellData.p_s = Vector(0.0, 0.0);
		// }
		//
		//
		//
		// cellData.u_s = cellData.p_s * (1 / reg(rho_s * cellData.eps_s));
		// /* //I don't thing this cause problem
		// if (cellData.eps_s > 0.7){
		// cellData.eps_s = 0.7; // sand packaging limit
		// }*/
		// cellData.p_g += (cellData.fluxP_g * TimeStep);
		//
		// cellData.rho_g += cellData.fluxRho_g * (TimeStep / reg(cellData.eps_g));
		//
		// cellData.eps_g = 1.0 - cellData.eps_s;
		//
		// cellData.u_g = cellData.p_g * (1 / reg(cellData.rho_g * cellData.eps_g));
		//
		// //cellData.T = cellData.T;
		//
		// cellData.p = cellData.rho_g * R_spec * T;
		//
		// }



		double MultiphaseFlow::FlowModulation(const Vertex<MeshType::meshDimension(), double>& x)
		@@ -338,6 +283,11 @@ void MultiphaseFlow::setupMeshData(const std::string& fileName){

		meshData.at(face).RightCellKoef = (face.getCenter() - rv).normEukleid() / dists.at(face);

		double sum = meshData.at(face).LeftCellKoef + meshData.at(face).RightCellKoef;

		meshData.at(face).LeftCellKoef /= sum;

		meshData.at(face).RightCellKoef /= sum;
		}


		@@ -431,16 +381,18 @@ void MultiphaseFlow::InitializePointData()
		}


		void MultiphaseFlow::ActualizePointData(const MeshDataContainer<FlowData, MeshType::meshDimension()>& data) {
		void MultiphaseFlow::UpdateVertexData(const MeshDataContainer<FlowData, MeshType::meshDimension()>& data) {


		// Initialize vector with zeros
		for(size_t i = 0; i < meshData.getDataByDim<0>().size(); i++) {
		meshData.getDataByDim<0>().at(i).u_g = {};
		meshData.getDataByDim<0>().at(i).u_s = {};
		}

		// can run in parallel without limitations
		for (const auto& vert : mesh.getVertices()) {

		#pragma omp for
		for (size_t vertIndex = 0; vertIndex < mesh.getVertices().size(); vertIndex++){
		const auto& vert = mesh.getVertices()[vertIndex];
		// initialize the vectors with zeros
		meshData[vert].u_g = {};
		meshData[vert].u_s = {};
		for (const auto& cellIndex : vertToCellCon.at(vert)) {

		const MeshType::Cell& cell = mesh.getCells().at(cellIndex);

Unstructured_mesh/multiphaseflow.h

+12 −6

Original line number	Diff line number	Diff line
		@@ -175,6 +175,8 @@ struct EdgeData {
		double fluxRho_g; //flux of "mass" over cell boundary
		};

		MAKE_ATTRIBUTE_TRAIT(EdgeData, fluxRho_g,fluxP_g,RightCellKoef,LeftCellKoef,n,Length,LengthOverDist)

		// Enumerator Type for expessing other cell and point position
		enum Type{
		INNER = 1,
		@@ -201,6 +203,7 @@ struct PointData
		double cellKoef;
		};

		MAKE_ATTRIBUTE_TRAIT(PointData, PointType, cellKoef, u_g, u_s)

		struct CellData{
		double invVolume;
		@@ -384,7 +387,7 @@ public:
		// // Computes velocities in the points on the start of the program
		void InitializePointData();
		//
		void ActualizePointData(const MeshDataContainer<FlowData, MeshType::meshDimension()> &data);
		void UpdateVertexData(const MeshDataContainer<FlowData, MeshType::meshDimension()> &data);
		//
		// // Calculates velocities during the computation step for each cell
		// void CalculateVertexData(Mesh::MeshCell &cell, const NumData<FlowData> &cellData);
		@@ -406,7 +409,6 @@ public:
		dataWriter.writeToStream(ofile, compData, writer);

		ofile.close();
		DBGMSG("Data eported");

		}
		};
		@@ -515,11 +517,15 @@ inline void MultiphaseFlow::ComputeFluxGas_inner(const FlowData &leftData, const
		Vector<ProblemDimension,double> fluxP_g = (-productOf_u_And_n) *
		(leftData.p_g * edgeData.LeftCellKoef + rightData.p_g * edgeData.RightCellKoef);




		// adding the element of pressure gradient
		fluxP_g -= (leftData.getPressure() * edgeData.LeftCellKoef +
		rightData.getPressure() * edgeData.RightCellKoef) * edgeData.n;



		fluxP_g *= edgeData.Length;


		@@ -556,7 +562,7 @@ inline void MultiphaseFlow::ComputeFluxGas_inflow(const FlowData& innerCellData,
		double productOf_u_And_n = (modulatedU * edgeData.n);

		/*
		** Preparing prtial derivatives wrt
		** Preparing partial derivatives wrt
		** edge orientation
		*/

		@@ -615,7 +621,7 @@ inline void MultiphaseFlow::ComputeFluxGas_outflow(const FlowData& innerCellData
		double eps_g_e = productOf_u_And_n > 0 ? innerCellData.getEps_g() : outFlow.getEps_g();

		/*
		** Preparing prtial derivatives wrt
		** Preparing partial derivatives wrt
		** edge orientation
		*/

		@@ -792,7 +798,7 @@ inline void MultiphaseFlow::ComputeFluxSolid_inflow(const FlowData& innerCellDat
		double productOf_u_And_n = (modulatedU * edgeData.n);

		/*
		** Preparing prtial derivatives wrt
		** Preparing partial derivatives wrt
		** edge orientation
		*/

		@@ -852,7 +858,7 @@ inline void MultiphaseFlow::ComputeFluxSolid_outflow(const FlowData& innerCellDa
		double eps_s_e = productOf_u_And_n > 0 ? innerCellData.eps_s : outFlow.eps_s;

		/*
		** Preparing prtial derivatives wrt
		** Preparing partial derivatives wrt
		** edge orientation
		*/
		// computing derivatives of velocity

src/UnstructuredMesh/MeshDataContainer/MeshDataIO/VTKMeshDataWriter.h

+2 −2

Original line number	Diff line number	Diff line
		@@ -147,7 +147,7 @@ class VTKMeshDataWriter {
		template<typename IndexType, typename Real>

		static void write(std::ostream& ost, const DataContainer<T, MeshDimension> &data, VTKMeshWriter<MeshDimension,IndexType, Real>& writer){
		DBGVAR(IsIndexable<typename DefaultIOTraits<T>::traitsType::template type<Index>>::value);
		//DBGVAR(IsIndexable<typename DefaultIOTraits<T>::traitsType::template type<Index>>::value);
		writeColumn<T, Index, IndexType, Real>(ost, data, writer);
		ost << std::endl;
		writeCellData<Traits<T, Types...>, Index + 1>::write(ost, data, writer);
		@@ -159,7 +159,7 @@ class VTKMeshDataWriter {
		struct writeCellData <Traits<T, Types...>, Index, std::enable_if_t<Index == Traits<T, Types...>::size() - 1>>{
		template< typename IndexType, typename Real>
		static void write(std::ostream& ost, const DataContainer<T, MeshDimension> &data, VTKMeshWriter<MeshDimension,IndexType, Real>& writer){
		DBGVAR(IsIndexable<typename DefaultIOTraits<T>::traitsType::template type<Index>>::value);
		//DBGVAR(IsIndexable<typename DefaultIOTraits<T>::traitsType::template type<Index>>::value);
		writeColumn<T, Index, IndexType, Real>(ost, data, writer);
		ost << std::endl;
		}