Fixing the particle writing method and providing access to fluid prop…

…erties.

FEM/rf_mfp_new.h

@@ -120,13 +120,13 @@ class CFluidProperties
 	double phi_0_tt(double T) const;
 	double EffectiveDiffusionCoef(int CIndex, double* variables = NULL); // AKS
 	void SetFemEleStd(FiniteElement::CFiniteElementStd* fem) { Fem_Ele_Std = fem; }
+	FiniteElement::CFiniteElementStd* Fem_Ele_Std;
 private:
 	int fluid_id; // specification of substance (NB JUN 09)
 	std::string name;
 	std::string cmpNm1, cmpNm2, cmpNm3, cmpNm4; // component name
 	int cmpN; // components number
 	// FEM
-	FiniteElement::CFiniteElementStd* Fem_Ele_Std;
 	long node; // OK4704
 	// Density
 	int density_model;

FEM/rf_random_walk.cpp

@@ -5980,8 +5980,9 @@ void DATWriteParticleFile(int current_time_step)
 	sprintf(now, "%i", current_time_step);
 	string nowstr = now;
 
-	string vtk_file_name =  pathJoin(defaultOutputPath, FileName);
-	vtk_file_name += "_RWPT_" + nowstr + ".particles.vtk";
+	string vtk_file_name = FileName + "_RWPT_";
+	vtk_file_name += nowstr;
+	vtk_file_name += ".particles.vtk";
 	fstream vtk_file (vtk_file_name.data(),ios::out);
 	vtk_file.setf(ios::scientific,ios::floatfield);
 	vtk_file.precision(12);
@@ -5990,21 +5991,19 @@ void DATWriteParticleFile(int current_time_step)
 	vtk_file.seekg(0L, ios::beg);
 
 	// Write Header
-	vtk_file << "# vtk DataFile Version 3.6.2"
-	         << "\n";
-	vtk_file << "Particle file: OpenGeoSys->Paraview. Current time (s) = " << RW->CurrentTime << "\n";
-	vtk_file << "ASCII"
-	         << "\n";
+	vtk_file << "# vtk DataFile Version 3.6.2" << "\n";
+	vtk_file << "Particle file: OpenGeoSys->Paraview. Current time (s) = " <<
+	RW->CurrentTime  << "\n";
+	vtk_file << "ASCII"  << "\n";
 	vtk_file << "\n";
-	vtk_file << "DATASET POLYDATA"
-	         << "\n"; // KR vtk_file << "DATASET PARTICLES"  << "\n";
-	vtk_file << "POINTS " << RW->numOfParticles << " double"
-	         << "\n";
+	vtk_file << "DATASET POLYDATA"  << "\n"; //KR vtk_file << "DATASET PARTICLES"  << "\n";
+	vtk_file << "POINTS " << RW->numOfParticles << " double" << "\n";
 
 	// Write particle locations
 	for (int i = 0; i < np; ++i)
-		vtk_file << RW->X[i].Now.x << " " << RW->X[i].Now.y << " " << RW->X[i].Now.z << "\n";
-
+	{
+		vtk_file << RW->X[i].Now.x << " " << RW->X[i].Now.y << " " << RW->X[i].Now.z <<	endl;
+	}
 	// KR add "vertices" block to create a correct VTK file
 	vtk_file << "VERTICES " << np << " " << (2 * np) << "\n";
 	for (int i = 0; i < np; ++i)
@@ -6013,23 +6012,21 @@ void DATWriteParticleFile(int current_time_step)
 	// Write particle identities
 	vtk_file << "\n";
 	vtk_file << "POINT_DATA " << RW->numOfParticles << "\n";
-	vtk_file << "SCALARS identity float 1"
-	         << "\n";
-	vtk_file << "LOOKUP_TABLE default"
-	         << "\n";
+	vtk_file << "SCALARS identity int 1" << "\n";
+	vtk_file << "LOOKUP_TABLE default" << "\n";
 	for (int i = 0; i < np; ++i)
 		vtk_file << RW->X[i].Now.identity << "\n";
 
 	// Write particle on_boundary or not
 	vtk_file << endl;
-	vtk_file << "SCALARS on_boundary float 1" << endl;
+    vtk_file << "SCALARS on_boundary int 1" << endl;
 	vtk_file << "LOOKUP_TABLE default" << endl;
 	for (int i = 0; i < np; ++i)
 		vtk_file << RW->X[i].Now.on_boundary << endl;
 
 	// Write particle vectors
 	vtk_file << endl;
-	vtk_file << "VECTORS velocity float" << endl;
+	vtk_file << "VECTORS velocity double" << endl;
 	for (int i = 0; i < np; ++i)
 		vtk_file << RW->X[i].Now.Vx << " " << RW->X[i].Now.Vy << " " << RW->X[i].Now.Vz << endl;
 

FEM/vtk.cpp

@@ -23,8 +23,9 @@
 #include "fem_ele_std.h" // for element velocity
 #include "makros.h"
 #include "rf_mmp_new.h"
-#include "FileTools.h"
 #include "OutputTools.h"
+#include "FileTools.h"
+#include "rf_random_walk.h"
 
 using namespace std;
 
@@ -1303,7 +1304,7 @@ bool CVTK::WriteElementValue(std::fstream& fin, bool output_data, COutput* out,
 			WriteDataArrayFooter(fin);
 	}
 
-	// MMP
+	//MMP VALUES OUTPUT
 	if (out->mmp_value_vector.size() > 0)
 	{
 		for (size_t i_mmp = 0; i_mmp < out->mmp_value_vector.size(); i_mmp++)
@@ -1329,13 +1330,7 @@ bool CVTK::WriteElementValue(std::fstream& fin, bool output_data, COutput* out,
 					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, false);
-						fem->Config();
-						fem->getShapeFunctionCentroid();
-						CMediumProperties* mmp = mmp_vector[ele->GetPatchIndex()];
-						double mat_value = ELEMENT_MMP_VALUES::getValue(mmp, mmp_id, i_e, 0, 1.0);
+                        double mat_value = getElementMMP(mmp_id, ele, m_pcs);
 						fin << mat_value << " ";
 					}
 					fin << "\n";
@@ -1365,8 +1360,10 @@ bool CVTK::WriteElementValue(std::fstream& fin, bool output_data, COutput* out,
 	// MFP
 	if (out->mfp_value_vector.size() > 0)
 	{
-		for (size_t i_mfp = 0; i_mfp < out->mfp_value_vector.size(); i_mfp++)
-		{
+  		  static double dbuff0[20];
+
+
+        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)
@@ -1385,16 +1382,18 @@ bool CVTK::WriteElementValue(std::fstream& fin, bool output_data, COutput* out,
 				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, false);
 						fem->Config();
-						fem->getShapeFunctionCentroid();
+                        fem->SetGaussPoint(0, gp_r, gp_s, gp_t);
+                        fem->ComputeShapefct(1, dbuff0);
 						CFluidProperties* mfp = mfp_vector[0];
-						mfp->SetFemEleStd(fem);
+                        mfp->Fem_Ele_Std = fem;
 						double mat_value = ELEMENT_MFP_VALUES::getValue(mfp, mfp_id);
 						fin << mat_value << " ";
 					}