Skip to content
Snippets Groups Projects
Commit 24e3b9a1 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Use CosseratVTKWriter also for Cosserat rods 

Throw out the AmiraMesh support instead.  I have not seen
anybody using Amira for years.
parent fadc99d3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/gfe/cosseratvtkwriter.hh
+ 28
− 0
View file @ 24e3b9a1
......@@ -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
}
}
......
dune/gfe/rodreader.hh deleted 100644 → 0
+ 0
− 119
View file @ fadc99d3
#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
dune/gfe/rodwriter.hh deleted 100644 → 0
+ 0
− 338
View file @ fadc99d3
#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
src/rod3d.cc
+ 3
− 2
View file @ 24e3b9a1
......@@ -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);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment