yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2378
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Wed 2010-07-21 10:00:51 +0200
message:
  1. Some VTK recorder fixes.
  2. cpm fixes (untested)
modified:
  core/Body.hpp
  doc/sphinx/introduction.rst
  doc/sphinx/user.rst
  pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp
  pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp
  pkg/dem/meta/ConcretePM.cpp
  pkg/dem/meta/ConcretePM.hpp


--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription

=== modified file 'core/Body.hpp'
--- core/Body.hpp	2010-06-08 21:22:07 +0000
+++ core/Body.hpp	2010-07-21 08:00:51 +0000
@@ -86,7 +86,8 @@
 		.def_readonly("id",&Body::id,"Unique id of this body") // overwrites automatic def_readwrite("id",...) earlier
 		.def_readonly("clumpId",&Body::clumpId,"Id of clump this body makes part of; invalid number if not part of clump; see :yref:`Body::isStandalone`, :yref:`Body::isClump`, :yref:`Body::isClumpMember` properties. \n\n This property is not meant to be modified directly from Python, use :yref:`O.bodies.appendClumped<BodyContainer.appendClumped>` instead.")
 		.def_readwrite("mat",&Body::material,"Shorthand for :yref:`Body::material`")
-		.add_property("isDynamic",&Body::isDynamic,"Deprecated synonym for :yref:`Body::dynamic`")
+		.add_property("isDynamic",&Body::isDynamic,&Body::setDynamic,"Deprecated synonym for :yref:`Body::dynamic` |ydeprecated|")
+		// .def_readwrite("dynamic",&Body::isDynamic,&Body::setDynamic,"Whether the body will be moved by forces.")
 		.def_readwrite("mask",&Body::groupMask,"Shorthand for :yref:`Body::groupMask`")
 		.add_property("isStandalone",&Body::isStandalone,"True if this body is neither clump, nor clump member; false otherwise.")
 		.add_property("isClumpMember",&Body::isClumpMember,"True if this body is clump member, false otherwise.")

=== modified file 'doc/sphinx/introduction.rst'
--- doc/sphinx/introduction.rst	2010-07-12 09:11:16 +0000
+++ doc/sphinx/introduction.rst	2010-07-21 08:00:51 +0000
@@ -371,6 +371,8 @@
 
 You will only rarely modify forces from Python; it is usually done in c++ code and relevant documentation can be found in the Programmer's manual.
 
+
+
 .. _function-components:
 
 Function components

=== modified file 'doc/sphinx/user.rst'
--- doc/sphinx/user.rst	2010-07-12 09:11:16 +0000
+++ doc/sphinx/user.rst	2010-07-21 08:00:51 +0000
@@ -1382,11 +1382,22 @@
 You can move between frames (snapshots that were saved) via the "Animation" menu. After setting the view angle, zoom etc to your satisfaction, the animation can be saved with *File/Save animation*.
 
 
+******************************
+Python specialties and tricks
+******************************
+
+.. perhaps turn this section into a list of FAQs on python as gathered from the yade-users list?
+
+
+ 
+
 **************
 Extending Yade
 **************
 
-new constitutive law
+* new particle shape
+* new constitutive law
+
 
 
 ****************

=== modified file 'pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp'
--- pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp	2010-07-08 14:58:21 +0000
+++ pkg/dem/Engine/GlobalEngine/VTKRecorder.cpp	2010-07-21 08:00:51 +0000
@@ -8,8 +8,8 @@
 #include<vtkUnstructuredGrid.h>
 #include<vtkXMLUnstructuredGridWriter.h>
 #include<vtkZLibDataCompressor.h>
-//#include<vtkXMLMultiBlockDataWriter.h>
-//#include<vtkMultiBlockDataSet.h>
+#include<vtkXMLMultiBlockDataWriter.h>
+#include<vtkMultiBlockDataSet.h>
 #include<vtkTriangle.h>
 #include<vtkLine.h>
 #include<yade/core/Scene.hpp>
@@ -96,13 +96,17 @@
 	spheresAngVelLen->SetNumberOfComponents(1);
 	spheresAngVelLen->SetName("angVelLen");		//Length (magnitude) of angular velocity
 	
-	vtkSmartPointer<vtkFloatArray> spheresStressVec = vtkSmartPointer<vtkFloatArray>::New();
-	spheresStressVec->SetNumberOfComponents(3);
-	spheresStressVec->SetName("stressVec");
+	vtkSmartPointer<vtkFloatArray> spheresNormalStressVec = vtkSmartPointer<vtkFloatArray>::New();
+	spheresNormalStressVec->SetNumberOfComponents(3);
+	spheresNormalStressVec->SetName("normalStress");
+
+	vtkSmartPointer<vtkFloatArray> spheresShearStressVec = vtkSmartPointer<vtkFloatArray>::New();
+	spheresShearStressVec->SetNumberOfComponents(3);
+	spheresShearStressVec->SetName("shearStress");
 	
-	vtkSmartPointer<vtkFloatArray> spheresStressLen = vtkSmartPointer<vtkFloatArray>::New();
-	spheresStressLen->SetNumberOfComponents(1);
-	spheresStressLen->SetName("stressLen");
+	vtkSmartPointer<vtkFloatArray> spheresNormalStressNorm = vtkSmartPointer<vtkFloatArray>::New();
+	spheresNormalStressNorm->SetNumberOfComponents(1);
+	spheresNormalStressNorm->SetName("normalStressNorm");
 	
 	vtkSmartPointer<vtkFloatArray> spheresMaterialId = vtkSmartPointer<vtkFloatArray>::New();
 	spheresMaterialId->SetNumberOfComponents(1);
@@ -182,18 +186,11 @@
 			line->GetPointIds()->SetId(0,I->getId1());
 			line->GetPointIds()->SetId(1,I->getId2());
 			intrCells->InsertNextCell(line);
-			if(recActive[REC_CPM]){		//For CPM model 
-				const CpmPhys* phys = YADE_CAST<CpmPhys*>(I->interactionPhysics.get());
-				intrForceN->InsertNextValue(phys->Fn);
-				float fs[3]={abs(phys->shearForce[0]),abs(phys->shearForce[1]),abs(phys->shearForce[2])};
-				intrAbsForceT->InsertNextTupleValue(fs);
-			} else {									//For all other models
-				const NormShearPhys* phys = YADE_CAST<NormShearPhys*>(I->interactionPhysics.get());
-				float fn[3]={abs(phys->normalForce[0]),abs(phys->normalForce[1]),abs(phys->normalForce[2])};
-				float fs[3]={abs(phys->shearForce[0]),abs(phys->shearForce[1]),abs(phys->shearForce[2])};
-				intrForceN->InsertNextTupleValue(fn);
-				intrAbsForceT->InsertNextTupleValue(fs);
-			}
+			const NormShearPhys* phys = YADE_CAST<NormShearPhys*>(I->interactionPhysics.get());
+			float fn[3]={abs(phys->normalForce[0]),abs(phys->normalForce[1]),abs(phys->normalForce[2])};
+			float fs[3]={abs(phys->shearForce[0]),abs(phys->shearForce[1]),abs(phys->shearForce[2])};
+			intrForceN->InsertNextTupleValue(fn);
+			intrAbsForceT->InsertNextTupleValue(fs);
 		}
 	}
 
@@ -233,10 +230,13 @@
 					spheresAngVelLen->InsertNextValue(angVel.norm());
 				}
 				if(recActive[REC_STRESS]){
-					const Vector3r& stress = bodyStates[b->getId()].normStress+bodyStates[b->getId()].shearStress;
-					float s[3] = { stress[0],stress[1],stress[2] };
-					spheresStressVec->InsertNextTupleValue(s);
-					spheresStressLen->InsertNextValue(stress.norm());
+					const Vector3r& stress = bodyStates[b->getId()].normStress;
+					const Vector3r& shear = bodyStates[b->getId()].shearStress;
+					float n[3] = { stress[0],stress[1],stress[2] };
+					float s[3] = { shear [0],shear [1],shear [2] };
+					spheresNormalStressVec->InsertNextTupleValue(n);
+					spheresShearStressVec->InsertNextTupleValue(s);
+					spheresNormalStressNorm->InsertNextValue(stress.norm());
 				}
 				
 				if (recActive[REC_CPM]){
@@ -294,8 +294,8 @@
 	vtkSmartPointer<vtkDataCompressor> compressor;
 	if(compress) compressor=vtkSmartPointer<vtkZLibDataCompressor>::New();
 
+	vtkSmartPointer<vtkUnstructuredGrid> spheresUg = vtkSmartPointer<vtkUnstructuredGrid>::New();
 	if (recActive[REC_SPHERES]){
-		vtkSmartPointer<vtkUnstructuredGrid> spheresUg = vtkSmartPointer<vtkUnstructuredGrid>::New();
 		spheresUg->SetPoints(spheresPos);
 		spheresUg->SetCells(VTK_VERTEX, spheresCells);
 		spheresUg->GetPointData()->AddArray(radii);
@@ -310,8 +310,9 @@
 			spheresUg->GetPointData()->AddArray(spheresAngVelLen);
 		}
 		if (recActive[REC_STRESS]){
-			spheresUg->GetPointData()->AddArray(spheresStressVec);
-			spheresUg->GetPointData()->AddArray(spheresStressLen);
+			spheresUg->GetPointData()->AddArray(spheresNormalStressVec);
+			spheresUg->GetPointData()->AddArray(spheresShearStressVec);
+			spheresUg->GetPointData()->AddArray(spheresNormalStressNorm);
 		}
 		if (recActive[REC_CPM]){
 			spheresUg->GetPointData()->AddArray(cpmDamage);
@@ -326,15 +327,17 @@
 
 		if (recActive[REC_MATERIALID]) spheresUg->GetPointData()->AddArray(spheresMaterialId);
 		
-		vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
-		if(compress) writer->SetCompressor(compressor);
-		string fn=fileName+"spheres."+lexical_cast<string>(scene->currentIteration)+".vtu";
-		writer->SetFileName(fn.c_str());
-		writer->SetInput(spheresUg);
-		writer->Write();
+		if(!multiblock){
+			vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
+			if(compress) writer->SetCompressor(compressor);
+			string fn=fileName+"spheres."+lexical_cast<string>(scene->currentIteration)+".vtu";
+			writer->SetFileName(fn.c_str());
+			writer->SetInput(spheresUg);
+			writer->Write();
+		}
 	}
+	vtkSmartPointer<vtkUnstructuredGrid> facetsUg = vtkSmartPointer<vtkUnstructuredGrid>::New();
 	if (recActive[REC_FACETS]){
-		vtkSmartPointer<vtkUnstructuredGrid> facetsUg = vtkSmartPointer<vtkUnstructuredGrid>::New();
 		facetsUg->SetPoints(facetsPos);
 		facetsUg->SetCells(VTK_TRIANGLE, facetsCells);
 		if (recActive[REC_COLORS]) facetsUg->GetCellData()->AddArray(facetsColors);
@@ -344,35 +347,41 @@
 		}
 		if (recActive[REC_MATERIALID]) facetsUg->GetCellData()->AddArray(facetsMaterialId);
 		if (recActive[REC_MASK]) facetsUg->GetCellData()->AddArray(facetsMask);
-		vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
-		if(compress) writer->SetCompressor(compressor);
-		string fn=fileName+"facets."+lexical_cast<string>(scene->currentIteration)+".vtu";
-		writer->SetFileName(fn.c_str());
-		writer->SetInput(facetsUg);
-		writer->Write();	
+		if(!multiblock){
+			vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
+			if(compress) writer->SetCompressor(compressor);
+			string fn=fileName+"facets."+lexical_cast<string>(scene->currentIteration)+".vtu";
+			writer->SetFileName(fn.c_str());
+			writer->SetInput(facetsUg);
+			writer->Write();	
+		}
 	}
+	vtkSmartPointer<vtkUnstructuredGrid> intrUg = vtkSmartPointer<vtkUnstructuredGrid>::New();
 	if (recActive[REC_INTR]){
-		vtkSmartPointer<vtkUnstructuredGrid> intrUg = vtkSmartPointer<vtkUnstructuredGrid>::New();
 		intrUg->SetPoints(intrBodyPos);
 		intrUg->SetCells(VTK_LINE, intrCells);
-		if (recActive[REC_CPM]){
-			 intrUg->GetCellData()->AddArray(intrForceN);
-			 intrUg->GetCellData()->AddArray(intrAbsForceT);
+		intrUg->GetCellData()->AddArray(intrForceN);
+		intrUg->GetCellData()->AddArray(intrAbsForceT);
+		if(!multiblock){
+			vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
+			if(compress) writer->SetCompressor(compressor);
+			string fn=fileName+"intrs."+lexical_cast<string>(scene->currentIteration)+".vtu";
+			writer->SetFileName(fn.c_str());
+			writer->SetInput(intrUg);
+			writer->Write();
 		}
-		vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
-		if(compress) writer->SetCompressor(compressor);
-		string fn=fileName+"intrs."+lexical_cast<string>(scene->currentIteration)+".vtu";
+	}
+
+	if(multiblock){
+		vtkSmartPointer<vtkMultiBlockDataSet> multiblockDataset = vtkSmartPointer<vtkMultiBlockDataSet>::New();
+		int i=0;
+		if(recActive[REC_SPHERES]) multiblockDataset->SetBlock(i++,spheresUg);
+		if(recActive[REC_FACETS]) multiblockDataset->SetBlock(i++,facetsUg);
+		if(recActive[REC_INTR]) multiblockDataset->SetBlock(i++,intrUg);
+		vtkSmartPointer<vtkXMLMultiBlockDataWriter> writer = vtkSmartPointer<vtkXMLMultiBlockDataWriter>::New();
+		string fn=fileName+lexical_cast<string>(scene->currentIteration)+".vtm";
 		writer->SetFileName(fn.c_str());
-		writer->SetInput(intrUg);
-		writer->Write();
+		writer->SetInput(multiblockDataset);
+		writer->Write();	
 	}
-
-	//vtkSmartPointer<vtkMultiBlockDataSet> multiblockDataset = vtkSmartPointer<vtkMultiBlockDataSet>::New();
-	//multiblockDataset->SetBlock(0, spheresUg );
-	//multiblockDataset->SetBlock(1, facetsUg );
-	//vtkSmartPointer<vtkXMLMultiBlockDataWriter> writer = vtkSmartPointer<vtkXMLMultiBlockDataWriter>::New();
-	//string fn=fileName+lexical_cast<string>(scene->currentIteration)+".vtm";
-	//writer->SetFileName(fn.c_str());
-	//writer->SetInput(multiblockDataset);
-	//writer->Write();	
 }

=== modified file 'pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp'
--- pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp	2010-07-08 14:58:21 +0000
+++ pkg/dem/Engine/GlobalEngine/VTKRecorder.hpp	2010-07-21 08:00:51 +0000
@@ -10,7 +10,8 @@
 		((bool,compress,false,"Compress output XML files [experimental]."))
 		((bool,skipFacetIntr,true,"Skip interactions with facets, when saving interactions"))
 		((bool,skipNondynamic,false,"Skip non-dynamic spheres (but not facets)."))
-		((string,fileName,"","Base file name; it will be appended with {spheres,intrs,facets}-243100.vtu depending on active recorders and step number (243100 in this case). It can contain slashes, but the directory must exist already."))
+		((bool,multiblock,false,"Use multi-block (``.vtm``) files to store data, rather than separate ``.vtu`` files."))
+		((string,fileName,"","Base file name; it will be appended with {spheres,intrs,facets}-243100.vtu (unless *multiblock* is ``True``) depending on active recorders and step number (243100 in this case). It can contain slashes, but the directory must exist already."))
 		((vector<string>,recorders,,"List of active recorders (as strings). Accepted recorders are:\n\n``all``\n\tSaves all possible parameters, except of specific. Default value.\n``spheres``\n\tSaves positions and radii (``radii``) of :yref:`spherical<Sphere>` particles.\n``id``\n\tSaves id's (field ``id``) of spheres; active only if ``spheres`` is active.\n``clumpId``\n\tSaves id's of clumps to which each sphere belongs (field ``clumpId``); active only if ``spheres`` is active.\n``colors``\n\tSaves colors of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``color``); only active if ``spheres`` or ``facets`` are activated.\n``mask``\n\tSaves groupMasks of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>` (field ``mask``); only active if ``spheres`` or ``facets`` are activated.\n``materialId``\n\tSaves materialID of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>`; only active if ``spheres`` or ``facets`` are activated.\n``velocity``\n\tSaves linear and angular velocities of spherical particles as Vector3 and length(fields ``linVelVec``, ``linVelLen`` and ``angVelVec``, ``angVelLen`` respectively``); only effective with ``spheres``.\n``facets``\n\tSave :yref:`facets<Facet>` positions (vertices).\n``stress``\n\tSaves stresses of :yref:`spheres<Sphere>` and of :yref:`facets<Facet>`  as Vector3 and length; only active if ``spheres`` or ``facets`` are activated.\n``cpm``\n\tSaves data pertaining to the :yref:`concrete model<Law2_Dem3DofGeom_CpmPhys_Cpm>`: ``cpmDamage`` (normalized residual strength averaged on particle), ``cpmSigma`` (stress on particle, normal components), ``cpmTau`` (shear components of stress on particle), ``cpmSigmaM`` (mean stress around particle); ``intr`` is activated automatically by ``cpm``.\n``intr``\n\tWhen ``cpm`` is used, it saves magnitude of normal (``forceN``) and shear (``absForceT``) forces.\n\n\tWithout ``cpm``, saves [TODO]\n``rpm``\n\tSaves data pertaining to the :yref:`rock particle model<Law2_Dem3DofGeom_RockPMPhys_Rpm>`: ``rpmSpecNum`` shows different pieces of separated stones, only ids.\n ``rpmSpecMass`` shows masses of separated stones.\n\n"))
 		((int,mask,0,"If mask defined, only bodies with corresponding groupMask will be exported. If 0 - all bodies will be exported.")),
 		/*ctor*/

=== modified file 'pkg/dem/meta/ConcretePM.cpp'
--- pkg/dem/meta/ConcretePM.cpp	2010-07-18 13:20:57 +0000
+++ pkg/dem/meta/ConcretePM.cpp	2010-07-21 08:00:51 +0000
@@ -245,24 +245,27 @@
 	// (1) getting intial equilibrium distance (not stored in ScGeom and distFactor not accessible from here)
 	// (2) applying contact forces back on particles
 	const Vector3r& pos1(Body::byId(I->getId1(),scene)->state->pos); const Vector3r& pos2(Body::byId(I->getId2(),scene)->state->pos);
+	Real refLength;
 	// just the first time
 	if(I->isFresh(scene)){
+		// done with real sphere radii
 		Real minRad=(contGeom->refR1<=0?contGeom->refR2:(contGeom->refR2<=0?contGeom->refR1:min(contGeom->refR1,contGeom->refR2)));
 		BC->crossSection=Mathr::PI*pow(minRad,2);
-		BC->refLength=(pos2-pos1).norm();
-		BC->kn=BC->crossSection*BC->E/BC->refLength;
-		BC->ks=BC->crossSection*BC->G/BC->refLength;
+		// scale sphere's radii to effective radii (intial equilibrium)
+		Real refLength=(pos2-pos1).norm(); Real distCurr=contGeom->radius1+contGeom->radius2;
+		contGeom->radius1*=refLength/distCurr; contGeom->radius2*=refLength/distCurr;
+		contGeom->penetrationDepth=0;
+		BC->kn=BC->crossSection*BC->E/refLength;
+		BC->ks=BC->crossSection*BC->G/refLength;
 	}
 	// shorthands
 	Real& epsN(BC->epsN);
 	Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& omegaThreshold(Law2_ScGeom_CpmPhys_Cpm::omegaThreshold); const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage); Real& omega(BC->omega); Real& sigmaN(BC->sigmaN);  Vector3r& sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs); // for python access
 	const bool& isCohesive(BC->isCohesive);
-	// FIXME: penetrationDepth does not account for interactionDetectionFactor!
-	epsN=(pos2-pos1).norm()/BC->refLength-1;
+	epsN=contGeom->penetrationDepth/refLength;
 	
-	// FIXME: sign?
 	epsT=contGeom->rotate(epsT);
-	epsT+=contGeom->shearIncrement()/BC->refLength; 
+	epsT+=contGeom->shearIncrement()/refLength; 
 
 	// simplified public model
 	epsN+=BC->isoPrestress/E;

=== modified file 'pkg/dem/meta/ConcretePM.hpp'
--- pkg/dem/meta/ConcretePM.hpp	2010-07-17 17:37:54 +0000
+++ pkg/dem/meta/ConcretePM.hpp	2010-07-21 08:00:51 +0000
@@ -138,7 +138,6 @@
 			((Real,epsPlSum,0,"cummulative shear plastic strain measure (scalar) on this contact"))
 			((bool,isCohesive,false,"if not cohesive, interaction is deleted when distance is greater than zero."))
 			((Vector3r,epsT,Vector3r::Zero(),"Total shear strain (either computed from increments with :yref:`ScGeom` or simple copied with :yref:`Dem3DofGeom`) |yupdate|"))
-			((Real,refLength,NaN,"Reference contact length (only used with :yref:`Law2_ScGeom_CpmPhys_Cpm`)"))
 			,
 			createIndex(); epsT=Fs=Vector3r::Zero(); Fn=0; omega=0;
 			,