diff --git a/FEM/CMakeLists.txt b/FEM/CMakeLists.txt index 50f664cc4..8b53dfcf0 100644 --- a/FEM/CMakeLists.txt +++ b/FEM/CMakeLists.txt @@ -48,6 +48,7 @@ set( HEADERS Stiff_Bulirsch-Stoer.h tools.h vtk.h + OutputTools.h ) set( SOURCES diff --git a/FEM/OutputTools.h b/FEM/OutputTools.h new file mode 100644 index 000000000..ad01848f3 --- /dev/null +++ b/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 + +#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 &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 &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_ diff --git a/FEM/fem_ele_std.h b/FEM/fem_ele_std.h index 0828cccc9..0e79b59d2 100644 --- a/FEM/fem_ele_std.h +++ b/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 diff --git a/FEM/rf_mfp_new.h b/FEM/rf_mfp_new.h index dbfe550fa..495831c27 100644 --- a/FEM/rf_mfp_new.h +++ b/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) diff --git a/FEM/vtk.cpp b/FEM/vtk.cpp index 27e6a939e..941a3ec1e 100644 --- a/FEM/vtk.cpp +++ b/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_mmpmmp_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(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_mfpmfp_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(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; }