diff --git a/dune/gfe/cosseratvtkwriter.hh b/dune/gfe/cosseratvtkwriter.hh
index 97239597566cf33e209608d489acf58253925f0d..4b61207d8cfa3e425e535fc79217d8194a74af99 100644
--- a/dune/gfe/cosseratvtkwriter.hh
+++ b/dune/gfe/cosseratvtkwriter.hh
@@ -173,6 +173,8 @@ public:
connectivitySize += ((vtkOrder==2) ? 6 : 3) * nE.second;
else if (nE.first.isHexahedron())
connectivitySize += ((vtkOrder==2) ? 20 : 8) * nE.second;
+ else if (nE.first.isLine())
+ connectivitySize += ((vtkOrder==2) ? 3 : 2) * nE.second;
else
DUNE_THROW(Dune::IOError, "Unsupported element type '" << nE.first << "' found!");
}
@@ -324,6 +326,32 @@ public:
connectivity[i++] = localView.index(5);
connectivity[i++] = localView.index(7);
connectivity[i++] = localView.index(6);
+#endif
+ }
+ }
+
+ if (element.type().isLine())
+ {
+ if (vtkOrder==2)
+ {
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+ connectivity[i++] = localIndexSet.index(0);
+ connectivity[i++] = localIndexSet.index(2);
+ connectivity[i++] = localIndexSet.index(1);
+#else
+ connectivity[i++] = localView.index(0);
+ connectivity[i++] = localView.index(2);
+ connectivity[i++] = localView.index(1);
+#endif
+ }
+ else // first order
+ {
+#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7)
+ connectivity[i++] = localIndexSet.index(0);
+ connectivity[i++] = localIndexSet.index(1);
+#else
+ connectivity[i++] = localView.index(0);
+ connectivity[i++] = localView.index(1);
#endif
}
}
diff --git a/dune/gfe/rodreader.hh b/dune/gfe/rodreader.hh
deleted file mode 100644
index 552df29e1899e2bbc0ca4a83e7ff6e5a229a1cf3..0000000000000000000000000000000000000000
--- a/dune/gfe/rodreader.hh
+++ /dev/null
@@ -1,119 +0,0 @@
-#ifndef ROD_READER_HH
-#define ROD_READER_HH
-
-#include <vector>
-
-#include <dune/common/exceptions.hh>
-#include <dune/common/fvector.hh>
-#include <dune/istl/bvector.hh>
-
-#include "rigidbodymotion.hh"
-
-class RodReader
-{
-public:
-
- /** \brief Read a spatial rod from file
- */
- static void readRod(std::vector<RigidBodyMotion<double,3> >& rod,
- const std::string& filename)
- {
-
- // /////////////////////////////////////////////////////
- // Load the AmiraMesh file
- // ////////////////////////////////////////////////
-
- AmiraMesh* am = AmiraMesh::read(filename.c_str());
-
- if(!am)
- DUNE_THROW(Dune::IOError, "Could not open AmiraMesh file: " << filename);
-
- int i, j;
-
- bool datafound=false;
-
- // get the translation
- AmiraMesh::Data* am_ValueData = am->findData("Nodes", HxFLOAT, 3, "Coordinates");
- if (am_ValueData) {
- datafound = true;
-
- if (rod.size()<am->nElements("Nodes"))
- DUNE_THROW(Dune::IOError, "When reading data from a surface field the "
- << "array you provide has to have at least the size of the surface!");
-
- float* am_values_float = (float*) am_ValueData->dataPtr();
-
- for (i=0; i<am->nElements("Nodes"); i++) {
- for (j=0; j<3; j++)
- rod[i].r[j] = am_values_float[i*3+j];
- }
-
- } else
- DUNE_THROW(Dune::IOError, "Couldn't find Coordinates: " << filename);
-
- // get the rotation
- Dune::BlockVector<Dune::FieldVector<double,3> > dir0(rod.size()),dir1(rod.size());
- am_ValueData = am->findData("Nodes", HxFLOAT, 3, "Directors0");
- if (am_ValueData) {
- datafound = true;
-
- if (rod.size()<am->nElements("Nodes"))
- DUNE_THROW(Dune::IOError, "When reading data from a surface field the "
- << "array you provide has to have at least the size of the surface!");
-
- float* am_values_float = (float*) am_ValueData->dataPtr();
-
- for (i=0; i<am->nElements("Nodes"); i++) {
- for (j=0; j<3; j++)
- dir0[i][j] = am_values_float[i*3+j];
- }
-
- } else
- DUNE_THROW(Dune::IOError, "Couldn't find Directors0: " << filename);
-
- am_ValueData = am->findData("Nodes", HxFLOAT, 3, "Directors1");
- if (am_ValueData) {
- datafound = true;
-
- if (rod.size()<am->nElements("Nodes"))
- DUNE_THROW(Dune::IOError, "When reading data from a surface field the "
- << "array you provide has to have at least the size of the surface!");
-
- float* am_values_float = (float*) am_ValueData->dataPtr();
-
- for (i=0; i<am->nElements("Nodes"); i++) {
- for (j=0; j<3; j++)
- dir1[i][j] = am_values_float[i*3+j];
- }
-
- } else
- DUNE_THROW(Dune::IOError, "Couldn't find Directors1: " << filename);
-
-
- for (i=0; i<rod.size(); i++) {
- dir0[i] /= dir0[i].two_norm();
- dir1[i] /= dir1[i].two_norm();
-
- // compute third director as the outer normal of the cross-section
- Dune::FieldVector<double,3> dir2;
- dir2[0] = dir0[i][1]*dir1[i][2]-dir0[i][2]*dir1[i][1];
- dir2[1] = dir0[i][2]*dir1[i][0]-dir0[i][0]*dir1[i][2];
- dir2[2] = dir0[i][0]*dir1[i][1]-dir0[i][1]*dir1[i][0];
- dir2 /= dir2.two_norm();
-
- // the rotation matrix corresponding to the directors
- Dune::FieldMatrix<double,3,3> rot(0);
- for (j=0;j<3;j++) {
- rot[j][0] = dir0[i][j];
- rot[j][1] = dir1[i][j];
- rot[j][2] = dir2[j];
- }
- rod[i].q.set(rot);
- rod[i].q.normalize();
- }
-
- std::cout << "Rod successfully read from: " << filename << std::endl;
-
- }
-};
-#endif
diff --git a/dune/gfe/rodwriter.hh b/dune/gfe/rodwriter.hh
deleted file mode 100644
index c8026ba23280757928df86b387250a0d59f33e54..0000000000000000000000000000000000000000
--- a/dune/gfe/rodwriter.hh
+++ /dev/null
@@ -1,338 +0,0 @@
-#ifndef ROD_WRITER_HH
-#define ROD_WRITER_HH
-
-#include <fstream>
-#include <vector>
-#include <ctime>
-
-#include <dune/common/exceptions.hh>
-#include <dune/common/fvector.hh>
-#include <dune/istl/bvector.hh>
-#include <dune/solvers/common/numproc.hh>
-
-#include "rigidbodymotion.hh"
-
-class RodWriter
-{
-public:
-
- static void writeBinary(const std::vector<RigidBodyMotion<double,3> >& rod,
- const std::string& filename)
- {
- FILE* fpRod = fopen(filename.c_str(), "wb");
- if (!fpRod)
- DUNE_THROW(SolverError, "Couldn't open file " << filename << " for writing");
-
- for (size_t j=0; j<rod.size(); j++) {
-
- for (int k=0; k<3; k++)
- fwrite(&rod[j].r[k], sizeof(double), 1, fpRod);
-
- for (int k=0; k<4; k++) // 3d hardwired here!
- fwrite(&rod[j].q[k], sizeof(double), 1, fpRod);
-
- }
-
- fclose(fpRod);
-
- }
-
-};
-
-/** \brief Write a planar rod
- */
-void writeRod(const std::vector<RigidBodyMotion<double,2> >& rod,
- const std::string& filename)
-{
- int nLines = rod.size() + 1 + 3*rod.size();
-
- // One point for each center line vertex and two for a little director at each vertex
- int nPoints = 3*rod.size();
-
- double directorLength = 1/((double)rod.size());
-
- // /////////////////////
- // Write header
- // /////////////////////
-
- time_t rawtime;
- time (&rawtime);
-
- std::ofstream outfile(filename.c_str());
-
- outfile << "# AmiraMesh 3D ASCII 2.0" << std::endl;
- outfile << std::endl;
- outfile << "# CreationDate: " << ctime(&rawtime) << std::endl;
- outfile << std::endl;
- outfile << std::endl;
- outfile << "define Lines " << nLines << std::endl;
- outfile << "nVertices " << nPoints << std::endl;
- outfile << std::endl;
- outfile << "Parameters {" << std::endl;
- outfile << " ContentType \"HxLineSet\"" << std::endl;
- outfile << "}" << std::endl;
- outfile << std::endl;
- outfile << "Lines { int LineIdx } @1" << std::endl;
- outfile << "Vertices { float[3] Coordinates } @2" << std::endl;
- outfile << std::endl;
- outfile << "# Data section follows" << std::endl;
- outfile << "@1" << std::endl;
-
- // ///////////////////////////////////////
- // write lines
- // ///////////////////////////////////////
-
- // The center axis
- for (size_t i=0; i<rod.size(); i++)
- outfile << i << std::endl;
-
- outfile << "-1" << std::endl;
-
- // The directors
- for (size_t i=0; i<rod.size(); i++) {
- outfile << rod.size()+2*i << std::endl;
- outfile << rod.size()+2*i+1 << std::endl;
- outfile << "-1" << std::endl;
- }
-
- // ///////////////////////////////////////
- // Write the vertices
- // ///////////////////////////////////////
-
- outfile << std::endl << "@2" << std::endl;
-
- // The center axis
- for (size_t i=0; i<rod.size(); i++)
- outfile << rod[i].r[0] << " " << rod[i].r[1] << " 0" << std::endl;
-
- // The directors
- for (size_t i=0; i<rod.size(); i++) {
-
- Dune::FieldVector<double, 2> director;
- director[0] = -cos(rod[i].q.angle_);
- director[1] = sin(rod[i].q.angle_);
- director *= directorLength;
-
- outfile << rod[i].r[0]+director[0] << " " << rod[i].r[1]+director[1] << " 0 " << std::endl;
- outfile << rod[i].r[0]-director[0] << " " << rod[i].r[1]-director[1] << " 0 " << std::endl;
-
- }
-
- std::cout << "Result written successfully to: " << filename << std::endl;
-
-}
-
-/** \brief Write a spatial rod
- */
-void writeRod(const std::vector<RigidBodyMotion<double,3> >& rod,
- const std::string& filename,
- double radius = 1.0)
-{
- int nPoints = rod.size();
-
- // /////////////////////
- // Write header
- // /////////////////////
-
- time_t rawtime;
- time (&rawtime);
-
- std::ofstream outfile(filename.c_str());
-
- outfile << "# AmiraMesh 3D ASCII 2.0" << std::endl;
- outfile << std::endl;
- outfile << "# CreationDate: " << ctime(&rawtime) << std::endl;
- outfile << std::endl;
- outfile << std::endl;
- outfile << "nNodes " << nPoints << std::endl;
- outfile << std::endl;
- outfile << "Parameters {" << std::endl;
- outfile << " ContentType \"Rod\"" << std::endl;
- outfile << "}" << std::endl;
- outfile << std::endl;
- outfile << "Nodes { float[3] Coordinates } @1" << std::endl;
- outfile << "Nodes { float[3] Directors0 } @2" << std::endl;
- outfile << "Nodes { float[3] Directors1 } @3" << std::endl;
- outfile << std::endl;
- outfile << "# Data section follows" << std::endl;
-
- // ///////////////////////////////////////
- // Write the center axis
- // ///////////////////////////////////////
- outfile << "@1" << std::endl;
-
- for (size_t i=0; i<rod.size(); i++)
- //outfile << rod[i].r[0] << " " << rod[i].r[1] << " " << rod[i].r[2] << std::endl;
- outfile << rod[i].r << std::endl;
-
-
- // ///////////////////////////////////////
- // Write the directors
- // ///////////////////////////////////////
- outfile << std::endl << "@2" << std::endl;
-
- for (size_t i=0; i<rod.size(); i++) {
- Dune::FieldVector<double,3> dir = rod[i].q.director(0);
- dir *= radius;
- outfile << dir << std::endl;
- }
-
- outfile << std::endl << "@3" << std::endl;
-
- for (size_t i=0; i<rod.size(); i++) {
- Dune::FieldVector<double,3> dir = rod[i].q.director(1);
- dir *= radius;
- outfile << dir << std::endl;
- }
-
- std::cout << "Result written successfully to: " << filename << std::endl;
-
-}
-
-/** \brief Write a spatial rod along with a strain or stress field
- */
-void writeRodElementData(Dune::BlockVector<Dune::FieldVector<double, 1> >& data,
- const std::string& filename)
-{
- // /////////////////////
- // Write header
- // /////////////////////
-
- time_t rawtime;
- time (&rawtime);
-
- std::ofstream outfile(filename.c_str());
-
- outfile << "# AmiraMesh 3D ASCII 2.0" << std::endl;
- outfile << std::endl;
- outfile << "# CreationDate: " << ctime(&rawtime) << std::endl;
- outfile << std::endl;
- outfile << std::endl;
- outfile << "define Segments " << data.size() << std::endl;
- outfile << std::endl;
- outfile << "Parameters {" << std::endl;
- outfile << " ContentType \"ScalarRodField\"" << std::endl;
- outfile << " Encoding \"OnSegments\"" << std::endl;
- outfile << "}" << std::endl;
- outfile << std::endl;
- outfile << "Segments { float Data } @1" << std::endl;
- outfile << std::endl;
-
- outfile << "# Data section follows" << std::endl;
- outfile << "@1" << std::endl;
-
- // ///////////////////////////////////////
- // write data
- // ///////////////////////////////////////
-
- for (size_t i=0; i<data.size(); i++)
- outfile << data[i] << std::endl;
-
- std::cout << "Result written successfully to: " << filename << std::endl;
-
-}
-
-void writeRodStress(Dune::BlockVector<Dune::FieldVector<double, 6> >& data,
- const std::string& basename)
-{
- Dune::BlockVector<Dune::FieldVector<double, 1> > scalarStress(data.size());
-
- // Extract separate stress values and write them
- for (int i=0; i<6; i++) {
-
- for (size_t j=0; j<data.size(); j++)
- scalarStress[j] = data[j][i];
-
- // create a name
- std::stringstream iAsAscii;
- iAsAscii << i;
-
- // write the scalar field
- writeRodElementData(scalarStress, basename + ".stress" + iAsAscii.str());
-
- }
-
-}
-
-/** \brief Write a spatial rod bundle
- */
-void writeRodBundle(const std::vector<std::vector<RigidBodyMotion<double,3> > >& rods,
- const std::string& filename,
- double radius = 1.0)
-{
- int nRods = rods.size();
-
- int nPoints(0);
- for (int i=0; i<nRods; i++)
- nPoints += rods[i].size();
-
- // /////////////////////
- // Write header
- // /////////////////////
-
- time_t rawtime;
- time (&rawtime);
-
- std::ofstream outfile(filename.c_str());
-
- outfile << "# AmiraMesh 3D ASCII 2.0" << std::endl;
- outfile << std::endl;
- outfile << "# CreationDate: " << ctime(&rawtime) << std::endl;
- outfile << std::endl;
- outfile << std::endl;
- outfile << "define Rods " << nRods << std::endl;
- outfile << "nNodes " << nPoints << std::endl;
- outfile << std::endl;
- outfile << "Parameters {" << std::endl;
- outfile << " ContentType \"RodBundle\"" << std::endl;
- outfile << "}" << std::endl;
- outfile << std::endl;
- outfile << "Rods { int Vertices } @1" << std::endl;
- outfile << "Nodes { float[3] Coordinates } @2" << std::endl;
- outfile << "Nodes { float[3] Directors0 } @3" << std::endl;
- outfile << "Nodes { float[3] Directors1 } @4" << std::endl;
- outfile << std::endl;
- outfile << "# Data section follows" << std::endl;
-
- // ///////////////////////////////////////
- // Write the center axis
- // ///////////////////////////////////////
- outfile << "@1" << std::endl;
-
- for (size_t i=0; i<rods.size(); i++)
- outfile << rods[i].size() << std::endl;
- //outfile << rod[i].r[0] << " " << rod[i].r[1] << " " << rod[i].r[2] << std::endl;
-
- outfile << std::endl << "@2" << std::endl;
- for (size_t i=0; i<rods.size(); i++)
- for (size_t j=0; j<rods[i].size(); j++)
- outfile << rods[i][j].r << std::endl;
-
-
- // ///////////////////////////////////////
- // Write the directors
- // ///////////////////////////////////////
- outfile << std::endl << "@3" << std::endl;
-
- for (size_t i=0; i<rods.size(); i++)
- for (size_t j=0; j<rods[i].size(); j++) {
- Dune::FieldVector<double,3> dir = rods[i][j].q.director(0);
- dir *= radius;
- outfile << dir << std::endl;
- }
-
- outfile << std::endl << "@4" << std::endl;
-
- for (size_t i=0; i<rods.size(); i++)
- for (size_t j=0; j<rods[i].size(); j++) {
- Dune::FieldVector<double,3> dir = rods[i][j].q.director(1);
- dir *= radius;
- outfile << dir << std::endl;
- }
-
- std::cout << "Result written successfully to: " << filename << std::endl;
-
-}
-
-
-#endif
diff --git a/src/rod3d.cc b/src/rod3d.cc
index 302a8d1e6dc3d0cc06a03bf2319853151bb1b4e8..6eda1ec68973059a28ca404d8f26def8842f32a2 100644
--- a/src/rod3d.cc
+++ b/src/rod3d.cc
@@ -12,9 +12,9 @@
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/norms/energynorm.hh>
+#include <dune/gfe/cosseratvtkwriter.hh>
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/geodesicdifference.hh>
-#include <dune/gfe/rodwriter.hh>
#include <dune/gfe/rotation.hh>
#include <dune/gfe/rodassembler.hh>
#include <dune/gfe/riemanniantrsolver.hh>
@@ -160,7 +160,8 @@ int main (int argc, char *argv[]) try
// //////////////////////////////
// Output result
// //////////////////////////////
- writeRod(x, resultPath + "rod3d.result");
+ CosseratVTKWriter<GridType>::write<FEBasis>(feBasis,x, resultPath + "rod3d-result");
+
BlockVector<FieldVector<double, 6> > strain(x.size()-1);
rodAssembler.getStrain(x,strain);