support MMP_VALUES and MFP_VALUES in PVD output

FEM/CMakeLists.txt

@@ -48,6 +48,7 @@ set( HEADERS
 	Stiff_Bulirsch-Stoer.h
 	tools.h
 	vtk.h
+	OutputTools.h
 )
 
 set( SOURCES

FEM/OutputTools.h

@@ -0,0 +1,148 @@
+/**
+ * \copyright
+ * Copyright (c) 2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef OUTPUTTOOLS_H_
+#define OUTPUTTOOLS_H_
+
+#include <iostream>
+
+#include "rf_mmp_new.h"
+#include "fem_ele_std.h"
+
+using namespace std;
+
+
+struct ELEMENT_MMP_VALUES
+{
+	static double getValue(CMediumProperties* mmp, int mmp_id, long i_e, double* gp, double theta)
+	{
+		double mat_value = .0;
+		switch (mmp_id)
+		{
+		case 0:
+			mat_value = mmp->Porosity(i_e, theta);
+			break;
+		case 1:
+			mat_value = mmp->PermeabilityTensor(i_e)[0];
+			break;
+		case 2:
+			mat_value = mmp->StorageFunction(i_e, gp, theta);
+			break;
+		default:
+			cout << "ELEMENT_MMP_VALUES::getValue(): no MMP values specified" << endl;
+			break;
+		}
+		return mat_value;
+	}
+
+	static int getMMPIndex(const std::string &mmp_name)
+	{
+		int mmp_id = -1;
+		if (mmp_name.compare("POROSITY") == 0)
+		{
+			mmp_id = 0;
+		}
+		else if (mmp_name.compare("PERMEABILITY") == 0)
+		{
+			mmp_id = 1;
+		}
+		else if (mmp_name.compare("STORAGE") == 0)
+		{
+			mmp_id = 2;
+		}
+		else
+		{
+			cout << "ELEMENT_MMP_VALUES::getMMPIndex(): no valid MMP values specified. " << mmp_name << endl;
+		}
+		return mmp_id;
+	}
+};
+
+struct ELEMENT_MFP_VALUES
+{
+	static double getValue(CFluidProperties* mfp, int mfp_id)
+	{
+		double mat_value = .0;
+		switch (mfp_id)
+		{
+		case 0:
+			mat_value = mfp->Density();
+			break;
+		case 1:
+			mat_value = mfp->Viscosity();
+			break;
+		default:
+			cout << "ELEMENT_MFP_VALUES: no MFP values specified" << endl;
+			break;
+		}
+		return mat_value;
+	}
+
+	static int getMFPIndex(const std::string &mfp_name)
+	{
+		int mfp_id = -1;
+		if (mfp_name.compare("DENSITY") == 0)
+		{
+			mfp_id = 0;
+		}
+		else if (mfp_name.compare("VISCOSITY") == 0)
+		{
+			mfp_id = 1;
+		}
+		else
+		{
+			cout << "ELEMENT_MFP_VALUES: no valid MFP values specified. " << mfp_name << endl;
+		}
+		return mfp_id;
+	}
+};
+
+inline double getElementMMP(int mmp_id, MeshLib::CElem* ele, CRFProcess* m_pcs)
+{
+	double gp[3] = { .0, .0, .0 };
+	double theta = 1.0;
+	int gp_r, gp_s, gp_t;
+	ele->SetOrder(false);
+	CFiniteElementStd* fem = m_pcs->GetAssember();
+	fem->ConfigElement(ele, m_pcs->m_num->ele_gauss_points, false);
+	fem->Config();
+	fem->SetGaussPoint(0, gp_r, gp_s, gp_t);
+	fem->ComputeShapefct(1);
+	CMediumProperties* mmp = mmp_vector[ele->GetPatchIndex()];
+	double val = ELEMENT_MMP_VALUES::getValue(mmp, mmp_id, ele->GetIndex(), gp, theta);
+	return val;
+}
+
+inline double getNodeMMP(int mmp_id, MeshLib::CFEMesh* m_msh, MeshLib::CNode* node, CRFProcess* m_pcs)
+{
+	const std::vector<size_t> &connected_ele_ids = node->getConnectedElementIDs();
+	double ele_avg = .0;
+	for (long i_e = 0; i_e < (long) connected_ele_ids.size(); i_e++)
+	{
+		MeshLib::CElem* ele = m_msh->ele_vector[connected_ele_ids[i_e]];
+		ele_avg += getElementMMP(mmp_id, ele, m_pcs);
+	}
+	ele_avg /= connected_ele_ids.size();
+	return ele_avg;
+}
+
+inline double getNodeElementValue(int ele_value_id, MeshLib::CFEMesh* m_msh, MeshLib::CNode* node, CRFProcess* m_pcs)
+{
+	const std::vector<size_t> &connected_ele_ids = node->getConnectedElementIDs();
+	double ele_avg = .0;
+	for (long i_e = 0; i_e < (long) connected_ele_ids.size(); i_e++)
+	{
+		MeshLib::CElem* ele = m_msh->ele_vector[connected_ele_ids[i_e]];
+		ele_avg += m_pcs->GetElementValue(ele->GetIndex(), ele_value_id);
+	}
+	ele_avg /= connected_ele_ids.size();
+	return ele_avg;
+}
+
+#endif // OUTPUTTOOLS_H_

FEM/fem_ele_std.h

@@ -286,7 +286,9 @@ class CFiniteElementStd : public CElement
 	// Primary as water head
 	bool HEAD_Flag;
 	//
+public:
 	void Config();
+protected:
 	//
 	double CalCoefMass();
 	// 25.2.2007 WW

FEM/rf_mfp_new.h

@@ -123,6 +123,7 @@ class CFluidProperties
 	double phi_r_dd (double rho, double T) const;
 	double phi_0_tt (double T) const;
 	double EffectiveDiffusionCoef(int CIndex, double* variables = NULL); //AKS
+	void SetFemEleStd(FiniteElement::CFiniteElementStd* fem) { Fem_Ele_Std = fem;}
 
 private:
 	int fluid_id;                         // specification of substance (NB JUN 09)

FEM/vtk.cpp

@@ -24,6 +24,7 @@
 #include "makros.h"
 #include "rf_mmp_new.h"
 #include "FileTools.h"
+#include "OutputTools.h"
 
 using namespace std;
 
@@ -1321,5 +1322,110 @@ bool CVTK::WriteElementValue(std::fstream &fin,
 		if (!useBinary || !output_data)
 			WriteDataArrayFooter(fin);
 	}
+
+	//MMP
+	if(out->mmp_value_vector.size() > 0)
+	{
+		for (size_t i_mmp=0; i_mmp<out->mmp_value_vector.size(); i_mmp++) {
+			const std::string &mmp_name = out->mmp_value_vector[i_mmp];
+			int mmp_id = ELEMENT_MMP_VALUES::getMMPIndex(mmp_name);
+			if (mmp_id<0) continue;
+
+			if (!useBinary || !output_data)
+				WriteDataArrayHeader(fin, this->type_Double, mmp_name, 0, str_format, offset);
+
+			if (output_data)
+			{
+				if (m_pcs==NULL) {
+					m_pcs = PCSGetFlow();
+				}
+
+				if (!this->useBinary)
+				{
+					fin << "          ";
+					for(long i_e = 0; i_e < (long)msh->ele_vector.size(); i_e++)
+					{
+						ele = msh->ele_vector[i_e];
+						double mat_value = getElementMMP(mmp_id, ele, m_pcs);
+						fin << mat_value << " ";
+					}
+					fin << "\n";
+				}
+				else
+				{
+					//OK411
+					write_value_binary<unsigned int>(fin, sizeof(int) * (long)msh->ele_vector.size());
+					for (long i_e = 0; i_e < (long)msh->ele_vector.size(); i_e++) {
+						ele = msh->ele_vector[i_e];
+						double mat_value = getElementMMP(mmp_id, ele, m_pcs);
+						write_value_binary(fin, mat_value);
+					}
+				}
+			} else {
+				offset += (long)msh->ele_vector.size() * sizeof(double) + SIZE_OF_BLOCK_LENGTH_TAG;
+			}
+
+			if (!useBinary || !output_data)
+				WriteDataArrayFooter(fin);
+		}
+	}
+
+	//MFP
+	if(out->mfp_value_vector.size() > 0)
+	{
+		for (size_t i_mfp=0; i_mfp<out->mfp_value_vector.size(); i_mfp++) {
+			const std::string &mfp_name = out->mfp_value_vector[i_mfp];
+			int mfp_id = ELEMENT_MFP_VALUES::getMFPIndex(mfp_name);
+			if (mfp_id<0) continue;
+
+			if (!useBinary || !output_data)
+				WriteDataArrayHeader(fin, this->type_Double, mfp_name, 0, str_format, offset);
+
+			if (output_data)
+			{
+				if (m_pcs==NULL) {
+					m_pcs = PCSGetFlow();
+				}
+
+				if (!this->useBinary)
+				{
+					fin << "          ";
+					int gp_r, gp_s, gp_t;
+					for(long i_e = 0; i_e < (long)msh->ele_vector.size(); i_e++)
+					{
+						ele = msh->ele_vector[i_e];
+						ele->SetOrder(false);
+						CFiniteElementStd* fem = m_pcs->GetAssember();
+						fem->ConfigElement(ele, m_pcs->m_num->ele_gauss_points, false);
+						fem->Config();
+						fem->SetGaussPoint(0, gp_r, gp_s, gp_t);
+						fem->ComputeShapefct(1);
+						CFluidProperties* mfp = mfp_vector[0];
+						mfp->SetFemEleStd(fem);
+						double mat_value = ELEMENT_MFP_VALUES::getValue(mfp, mfp_id);
+						fin << mat_value << " ";
+					}
+					fin << "\n";
+				}
+				else
+				{
+					//OK411
+					write_value_binary<unsigned int>(fin, sizeof(int) * (long)msh->ele_vector.size());
+					for (long i_e = 0; i_e < (long)msh->ele_vector.size(); i_e++) {
+						ele = msh->ele_vector[i_e];
+						CFluidProperties* mfp = mfp_vector[0];
+						double mat_value = ELEMENT_MFP_VALUES::getValue(mfp, mfp_id);
+						write_value_binary(fin, mat_value);
+					}
+				}
+			} else {
+				offset += (long)msh->ele_vector.size() * sizeof(double) + SIZE_OF_BLOCK_LENGTH_TAG;
+			}
+
+			if (!useBinary || !output_data)
+				WriteDataArrayFooter(fin);
+		}
+	}
+
 	return true;
 }