[Date Prev][Date Next][Thread Prev][Thread Next]
[Date Index]
[Thread Index]
- From: christof@...
- Date: Thu, 26 Oct 2006 17:41:58 +0200
Sourcen.
In section.cpp so ab Zeile 225.
#include "Section.h"
#include "Archive.h"
#include "vtbutil.h"
#include <math.h>
#pragma warning(disable:4996)
const int Section::version = 154;
void Element::Print(int r)
{
char buf[255];
sprintf(buf, " %2d %6.2f %6.3f %6.2f\n", r+1, b, t, a);
cout << buf;
}
void Section::RW(vtbArchive& a)
{
a.Line(qw_id);
if (_loadedVersion == 153)
a.InfoLine("*** Typ ns nk nk1 nk2 nk3 Kraft Laenge nla nlas nlak nps ndb FG nz nzl nzb nzd ***");
else
a.InfoLine("*** Typ ns nk nk1 nk2 nk3 Kraft Laenge nla nlas nlak nps ndb nsf FG nz nzl nzb nzd ***");
a.rw(itype);
a.rw(ns);
a.rw(nk);
a.rw(nk1);
a.rw(nk2);
a.rw(nk3);
a.rw(forceUnit);
a.rw(lengthUnit);
a.rw(nla);
a.rw(nlas);
a.rw(nlak);
a.rw(nps);
a.rw(ndb);
if (_loadedVersion == 154)
a.rw(nsf);
a.rw(dof.letters);
a.rw(sv.nz);
a.rw(sv.nzl);
a.rw(sv.nzb);
a.rw(sv.nzd);
a.StopLine();
a.InfoLine("*** E-Modul G-MODUL MUE ***");
a.rw(eMod);
a.rw(gMod);
a.rw(nue);
a.StopLine();
a.InfoLine("*** Querschnittserlaeuterung ***");
a.rw(description);
a.InfoLine("*** Anfangs- Endknoten ***");
a.rw(ns, iak);
a.rw(ns, iek);
a.InfoLine("*** Knotenkoordinaten, Scheiben ***");
a.rw(nk, ky);
a.rw(nk, kz);
a.rw(nk, nsch);
a.InfoLine("*** Scheiben b_r t_r alpha_r K_r D_r ***");
a.rw(ns, br);
a.rw(ns, tr);
a.rw(ns, al);
a.rw(ns, kr);
// a.InfoLine("*** Scheiben ***");
// a.rw(element);
// a.InfoLine("*** Knoten ***");
// a.rw(node);
if (nk1 > 0){
a.InfoLine("*** Lagerung der Randknoten ***");
a.rw(nk1,randknoten_nr);
a.rw(nk1,randlager_typ);
a.rw(nk1,randlager_winkel);
}
else {
a.InfoLine("*** Keine Randlagerungen ***");
}
if (nps > 0){
a.InfoLine("*** Pendelstaebe: Knoten, Winkel ***");
a.rw(nps,ips);
a.rw(nps,al_ps);
a.rw(nps,c_ps);
a.StopLine();
}
else {
a.InfoLine("*** Keine Pendelstaebe an Innenknoten ***");
}
if (ndb > 0){
a.InfoLine("*** Drehfedern: Scheibe, Steifigkeit ***");
a.rw(ndb,idb);
a.rw(ndb,c_theta);
a.rw(ndb,c_ps);
a.StopLine();
}
else {
a.InfoLine("*** Keine Drehfedern ***");
}
if (_loadedVersion == 154)
{
if (nsf > 0){
a.InfoLine("*** Schubfedern: Knoten, Winkel, Steifigkeit ***");
a.rw(nsf,schubfeder_r);
a.rw(nsf,schubfeder_alpha);
a.rw(nsf,schubfeder_c);
a.StopLine();
}
else {
a.InfoLine("*** Keine Schubfedern ***");
}
}
if (nla > 0)
{
a.InfoLine("*** Lagerelemente ***");
if (nlak > 0){
a.InfoLine("*** Knoten Nr. Lagertyp Winkel Steifigkeit ***");
a.rw(nlak,knlag_nr);
a.rw(nlak,knlag_typ);
a.rw(nlak,knlag_al);
a.rw(nlak,knlag_c);
a.StopLine();
}
if (nlas > 0){
a.InfoLine("*** Scheibennr. Steifigkeit ***");
a.rw(nlas,slag_nr);
a.rw(nlas,slag_c);
a.StopLine();
}
}
a.InfoLine("*** Schwerpunkt, Hauptachsen ***");
a.rw(sv.centerOfGravity.x);
a.rw(sv.centerOfGravity.y);
a.rw(sv.mainAxisAngle[0]);
a.rw(sv.mainAxisAngle[1]);
a.StopLine();
a.InfoLine("*** Zustandsattribute ***");
a.rw(sv.nz, sv.modeAttributeEnum);
a.StopLine();
}
void SectionValues::RW(vtbArchive& a)
{
a.rw(qwReports);
a.rw(U, "*** U ***");
a.rw(F_s, "*** F_s ***");
a.rw(F_sq, "*** F_sq ***");
a.rw(F_th, "*** F_th ***");
a.rw(V, "*** V ***");
a.rw(W, "*** W ***");
a.rw(Sch_fl_a, "*** Schubfluss a ***");
a.rw(Sch_fl_e, "*** Schubfluss e ***");
a.rw(Sch_kr, "*** Schubkraefte ***");
a.rw(Phi, "*** Phi ***");
DegreeOfFreedom dof;
if (dof.parabolicMoments)
{
a.rw(M_p, "*** M_p ***");
}
if (dof.shearDeformation)
{
a.rw(Fs_q, "*** Fs_q ***");
a.rw(Fs_d, "*** Fs_d ***");
a.rw(Sch_kr_ga, "*** Schubkraefte aus gamma ***");
}
a.rw(Ck_mem, "*** C_mem ***");
a.rw(Ck, "*** C ***");
a.rw(D, "*** D ***");
a.rw(D1, "*** D1 ***");
a.rw(D2, "*** D2 ***");
//if (nsf)
{
a.rw(D3, "*** D3 ***");
}
a.rw(B, "*** B ***");
}
#include <sstream>
#include <iostream>
void Section::Load(std::string& filename)
{
vtbArchive a(filename, vtbArchive::READ);
std::string l = a.GetLine();
if (l.substr(0,4) != "****")
return;
_loadedVersion = int (Val1<float>(l.substr(32, 4)) * 100.0f + 0.5f);
RW(a);
}
void Section::Save(std::string& filename)
{
vtbArchive a(filename, vtbArchive::WRITE);
char buf[128];
sprintf(buf, "**** VTB-Datei, Formatversion = %4.2f ****", 0.01f * version);
std::string b(buf);
a.Line(b);
_loadedVersion = version;
RW(a);
}
void Section::SearchElements2(int iNode, int& i1, int& i2)
{
// Searches the elemnents, which are linked to node 'iNNode' and
// put their indices in i1 an i2.
const int EINS = 1; // Because of OPTION BASE 1 in GFA/QED
i1=-1, i2=-1;
for (int r=0; r<ns; ++r){
if (iak.at(r)-EINS == iNode)
i1 = r;
if (iek.at(r)-EINS == iNode)
i2 = r;
}
}
bool Section::IsCorner(int iNode)
{
// Gives true, if the node iNode is a 'Corner'. Which means:
// 1. must exactly have two plates linked
// 2. must have an difference-angle <> 0
const double delMin = 10e-3;
if (nsch.at(iNode) == 2) //
{
int i1, i2;
SearchElements2(iNode, i1, i2);
double deltaAlpha = al.at(i2) - al.at(i1);
return bool(abs(deltaAlpha) > delMin);
}
else
return false;
}
void Section::MakeRoundedCorners(double radius, int nAdded, bool convex, vector<int> divElements)
{
// Changes the section in some ways:
// - Corners are rounded by inserting small elemnts (on a circle-path)
// variables: nAdded - how many elements to insert
// convex - curvature in which direction
// radius - of the assumed circle-path
// - elements are subdivided (optional). For each original element there
// is a number in divElements, which will - if<>1 - be taken to divide the element.
// Strategy: walk down the original elements br[] and al[] and create
// a set of new elements bb and aa.
// Afterwards
if (itype == OPEN)
{
// New arrays for generated elements:
vector<double> aa; // angle
vector<double> bb; // width
for (int r=1; r<nk; ++r)// walk down the inner nodes
{
if (!IsCorner(r))
{
int& dv = divElements[r-1];
if (dv != 1)
{
double divBr = (br[r-1])/dv;
for (int j=0; j<dv; ++j) {
aa.push_back(al[r-1]);
bb.push_back(divBr);
}
}
else {
aa.push_back(al[r-1]);
bb.push_back(br[r-1]);
}
}
else
{
RoundedCornerData d;
CalcRoundedCorner(radius, nAdded, r, d);
// (1) add the element before the corner
int& dv = divElements[r-1];
if (dv != 1)
{ // if divided
double divBr = (br[r-1] - d.shorteningB)/dv;
for (int j=0; j<dv; ++j) {
aa.push_back(al[r-1]);
bb.push_back(divBr);
}
}
else
{ // undevided, just addd the shortened element
aa.push_back(al[r-1]);
bb.push_back(br[r-1] - d.shorteningB);
}
// (2) add the elements created for the corner
double a, da;
if (convex){
a = al[r-1];
da = -d.deltaAlpha;
} else {
a = al[r];
da = d.deltaAlpha;
}
a += da/2.0;
for (int i=0; i<d.nAddedElements; ++i){
aa.push_back(a);
bb.push_back(d.b);
a += da;
}
// (3) the element AFTER the corner is shortend:
// (it is added in the next step)
br[r] -= d.shorteningB;
}
}
// Now take the new vectors into the section-data
al = aa;
br = bb;
// ... and calculate the rest (the depenant values)
CalcOpenSection();
}
else
{
// General section (branched etc.)
// This has to be done..............................
int oldnk = nk;
for (int r=0; r<oldnk; ++r) {
if (IsCorner(r)) {
MakeRoundedCorner(radius, r);
}
}
}
// Beispiel für Winkel :
// 90-Grad-Ecke. 3 Scheiben. al0=90, da=30 da/2=15
// al = 0 15 45 75 90 Elements-Alpha
// 15 30 30 15 Delta-Alpha
}
void Section::CalcRoundedCorner(double radius, int nAdded, int iNode, RoundedCornerData& rcd)
{
// Calculates the characteristical data for a rounded corner
// and returns it in the element 'rcd'.
rcd.nAddedElements = nAdded;
int i1, i2;
SearchElements2(iNode, i1, i2);
double da = al[i2] - al[i1]; // Difference-Angle of main elements
rcd.deltaAlpha = da/rcd.nAddedElements; // Angel-Increment of inserted elements
rcd.shorteningB = abs(radius*tand(da/2.0)); // shorten the attached elements
if (rcd.shorteningB > br[i1]/2 || rcd.shorteningB > br[i2]/2 ) {
// error("Elements to short for radius") // todo... errormessage
return;
}
rcd.b = abs(2.0*radius*sind(rcd.deltaAlpha/2.0)); // width of new elements
}
void Section::MakeRoundedCorner(double radius, int iNode)
{
int n = 4;
int i1, i2;
SearchElements2(iNode, i1, i2);
double deltaAlpha = al[i1] - al[i2]; // Angle of main elements
double incAlpha = deltaAlpha/n; // Angel-Increment of inserted elements
if (deltaAlpha < 0)
deltaAlpha += 360.0;
double bCutoff = abs(radius*tand(deltaAlpha/2.0)); // shorten the attached elements
if (bCutoff > br[i1]/2 || bCutoff > br[i2]/2 ) {
// error("Elements to short for radius")
return;
}
double s = abs(2.0*radius*sind(incAlpha/2.0)); // width of new elements
}
void Section::AddElement(double b , double a)
{
al.push_back(a);
br.push_back(b);
}
void Section::CalcOpenSection()
{
// The topology of a simple section (open, unbranched) is completely defined
// by the width and angles of the elements.
// Here this values are taken and all other dependant values are calculated
// from that.
const int EINS = 1; // because GFA/QED-Files start indexing with 1, this is the correction
ns = br.size();
nk = ns + 1;
nk1 = 2;
nk2 = ns-1;
nk3 = 0;
iak.clear();
iek.clear();
for (int i=0; i<ns; ++i) {
iak.push_back(i + EINS);
iek.push_back(i+1 + EINS);
}
nsch.assign(nk,2);
nsch.front() = 1;
nsch.back() = 1;
ky.clear();
kz.clear();
ky.push_back(0.0);
kz.push_back(0.0);
for (int i=0; i<ns; ++i)
{
double y = ky.back();
double z = kz.back();
y -= cosd(al.at(i))*br.at(i);
z -= sind(al.at(i))*br.at(i);
ky.push_back(y);
kz.push_back(z);
}
// int ymin=
for (int i=0; i<nk; ++i)
{
//.......
}
randknoten_nr.clear();
randlager_winkel.clear();
randlager_typ.clear();
randknoten_nr.push_back(0 + EINS);
randknoten_nr.push_back(nk-1 + EINS);
randlager_typ.push_back(0);
randlager_typ.push_back(0);
randlager_winkel.push_back(0.0);
randlager_winkel.push_back(0.0);
}
/*
Weitermachen:
*/
Section::Section()
{
ns = 0;
nk = 0;
eMod = 21000.0;
nue = 0.3;
gMod = eMod/(2.0*(1.0+nue));
itype = OPEN;
nla = 0;
nlas = 0;
nlak = 0;
nsf = 0;
ndb = 0;
nps = 0;
forceUnit = "kN";
lengthUnit = "mm";
dof.letters = "u";
}
void Section::SetT(double t)
{
tr.assign(ns, t);
const double k = eMod*pow(t,3)/(12.0*(1-nue*nue));
kr.assign(ns, k);
}
SectionValues::SectionValues()
{
nz = 0;
centerOfGravity.x = 0.0;
centerOfGravity.y = 0.0;
mainAxisAngle[0] = 0.0;
mainAxisAngle[1] = 0.0;
}
void Section::Underline(string a, char c, bool above)
{
if (above)
cout << string(a.size(),c) << endl;
cout << a << endl;
cout << string(a.size(),c) << endl;
}
void Section::Print()
{
TransformToNewFormat();
Underline("Querschnittsbeschreibung",'=',true);
for (int i=0; i<description.size(); ++i)
cout << description[i] << endl;
cout << "Der Querschnitt besitzt "<<nk<<" Knoten und "<<ns<<" Scheiben." << endl;
if(IsScheibenzug()) {
if (IsClosed())
cout << "Es ist ein geschlossener Scheibenzug" << endl;
else
cout << "Es ist ein offener Scheibenzug" << endl;
}
else
cout << "Es gibt "<<nk3<<" Verzweigungsknoten und "<<nk1<<" Endknoten." << endl;
cout << endl;
cout << "Einheiten: " << forceUnit << " und " << lengthUnit << endl;
cout << endl;
if(IsScheibenzug()) {
cout << "Scheiben:" ;
cout << " r b_r t_r alpha_r " << endl;
for (int r=0; r<ns; ++r)
element[r].Print(r);
}
cout << endl;
}
void Section::TransformToNewFormat()
{
element.reserve(ns);
for (int i=0; i<ns; ++i)
{
Element e;
e.a = al[i];
e.b = br[i];
e.k = kr[i];
e.t = tr[i];
e.ia = iak[i];
e.ie = iek[i];
element.push_back(e);
}
node.reserve(ns);
for (int i=0; i<nk; ++i)
{
Node n;
n.y = ky[i];
n.z = kz[i];
n.nElements = nsch[i];
node.push_back(n);
}
}
#include <vector>
#include "MatVec.h"
//#include "Archive.h"
using namespace std;
class vtbArchive;
//namespace gbt {
/************************************************************************
* Structures describing
************************************************************************/
struct Element
{
double b; // width
double t; // thickness
double a; // angle, measured counterclockwise from three o'clock.
// unit: degree
double k; // plate-stiffness, calculated from e-modulus and t.
int ia; // Index of first node
int ie; // and second (ia->ie gives the orientiation for
// the local coordinate s)
void Print(int r);
};
// Coordinate-pair
struct DPoint
{
double x,y;
};
struct Node
{
double y,z;
// oder: DPoint k;
int nElements; // number of linked elements
};
/************************************************************************
* Structures describing constraints and boundary conditions
************************************************************************/
// for boundary nodes
struct BoundaryConstraint
{
int type; // type of constraint
int i; // index of node
double al; // angle
// These constraints are yet limited to rigid supports
// double c; // stiffness of spring ( 0 means
};
// for inner nodes
struct NodeConstraint
{
// int type; // type of constraint: always: translational spring
int i; // index of node
double al; // angle
double c; // stiffness of spring ( 0 means
};
// for elements
struct RotationalSpring // ElementConstraint
{
int i; // index of element
double c; // stiffness
};
// general constraint, not yet used
struct NodeSpring
{
int type; // type of constraint
int i; // index of node
double al; // angle
double c; // stiffness of spring
};
/************************************************************************
* Degrees of freedom, on which the section and -values are based
************************************************************************/
struct DegreeOfFreedom
{
std::string letters; // DOF, expressed in a string of letters
bool warping;
bool parabolicMoments;
bool shearDeformation;
};
// Todo.....:
enum SectionType {OPEN=1, CLOSED, I_SIMPLE, I_SYMM, I_GENERAL, BRANCHED, PLATE};
// Characterisation of the modes
struct ModeAttribute
{
int type; // 1=elongation 2=bending 3=torsion 4=distorsion
int symmetry; // 0=none 1=symm. 2=antisymm.
};
struct SectionValues;
/************************************************************************
* Section
************************************************************************/
struct SectionValues
{
SectionValues();
//-----------------------------------------------------------------------
// General information
//-----------------------------------------------------------------------
string qwId; // Identification of the generating program
vector<ModeAttribute> modeAttribute; // Mode-Characterisation
vector<int> modeAttributeEnum; // Mode-Characterisation
DPoint centerOfGravity; //
double mainAxisAngle[2]; //
string qwReports; // reports coming from the QW-program
//-----------------------------------------------------------------------
// Mode numbers
//-----------------------------------------------------------------------
int nz; // number of orthogonal modes
int nzl; // number of longitudinal modes (0 or 1)
int nzb; // number of bending modes (0,1 or 2)
int nzd; // number of rotaional modes (0 or 1)
// nz - (nzb+nzl+nzd) are the distortional modes
//-----------------------------------------------------------------------
// Modal deformations and stresses
//-----------------------------------------------------------------------
DMatrix U; // warping ordinates
DMatrix F_s; // element displ. in s-direction
DMatrix F_sq; // ..and s-bar-direction
DMatrix F_th; // element rotation (counterclockwise!!!)
DMatrix V, W; // Displacement of nodes
DMatrix Sch_fl_a; // shearflow
DMatrix Sch_fl_e; //
DMatrix Sch_kr; //
DMatrix Phi; //
// If option 'parabolic moments m_s'
DMatrix M_p; // parabolic transverse-moments
// If option 'shear & eps_s'
DMatrix Fs_q, Fs_d;
DMatrix Sch_kr_ga;
//-----------------------------------------------------------------------
// Stiffness parameters
//-----------------------------------------------------------------------
DVector Ck; // warping-resistance
DVector Ck_mem; // warping-resistance, only membrane part
DMatrix D; // D-matrix (
DMatrix D1; // D-matrix (part 1)
DMatrix D2; // D-matrix (part 2)
DMatrix D3; // D-matrix (part 3) ???????????????????????????????????
DMatrix B; // B-matrix (transverse bending stiffness)
//-----------------------------------------------------------------------
//
//-----------------------------------------------------------------------
void RW(vtbArchive& a);
};
struct Section
{
Section();
void Load(std::string& filename);
void Save(std::string& filename);
int _loadedVersion;
const static int version;
std::string qw_id;
SectionType type; // topologic type of section
int itype; // topologic type of section
//-----------------------------------------------------------------------
// Geometry, topology
//-----------------------------------------------------------------------
// number of ...
int ns; // elements
int nk; // nodes
int nk1; // nodes linked to one single element (boundary-node)
int nk2; // nodes linked to two elements (normal inner-nodes)
int nk3; // nodes linked to three elements (branch-nodes)
// vector<DPoint> node; // Vector of nodes
vector<Node> node; // Vector of nodes
vector<Element> element; // and elements
vector<int> iak;
vector<int> iek;
vector<double> ky;
vector<double> kz;
vector<int> nsch;
vector<double> br;
vector<double> tr;
vector<double> al;
vector<double> kr;
//-----------------------------------------------------------------------
// Material
//-----------------------------------------------------------------------
double eMod; // Youngs modulus
double nue; // Poisson-const.
double gMod; // Shear modulus (calculated from eMod and nue).
// e g my
//-----------------------------------------------------------------------
// Constraints
//-----------------------------------------------------------------------
int nps; // number of node-supports (at intermediate nodes)
int ndb; // number of rotational springs
int nsf; // number of shear-springs
vector<BoundaryConstraint> boundarySupport; // "randlager_", size is nk1
vector<NodeSpring> nodeSupport; // "ps"
vector<RotationalSpring> rotSpring; // "db"
int nlas; // number of nodal constraints
int nlak; // number of elment constraints
int nla; // sum of constraints (upper two)
vector<int> randknoten_nr;
vector<int> randlager_typ;
vector<double> randlager_winkel;
vector<int> ips;
vector<double> al_ps;
vector<double> c_ps;
vector<int> idb;
vector<double> c_theta;
vector<int> schubfeder_r;
vector<double> schubfeder_alpha;
vector<double> schubfeder_c;
vector<int> knlag_nr;
vector<int> knlag_typ;
vector<double> knlag_al;
vector<double> knlag_c;
vector<int> slag_nr;
vector<double> slag_c;
//-----------------------------------------------------------------------
// Additional information
//-----------------------------------------------------------------------
string forceUnit; // units
string lengthUnit;
vector<string> description; // arbitrary comments by the creator of the section
DegreeOfFreedom dof;
SectionValues sv;
void RW(vtbArchive& a);
//----------------------------------------
struct RoundedCornerData
{
int nAddedElements;
double shorteningB;
double b; // width of added elements
double deltaAlpha; // difference Angle betwwen added elements (half angle at the borders)
};
bool IsCorner(int i);
void SearchElements2(int iNode, int& i1, int& i2);
void MakeRoundedCorners(double radius, int nAdded, bool convex, vector<int> divElements);
void MakeRoundedCorner(double radius, int iNode);
void CalcRoundedCorner(double radius, int nAdded, int iNode, RoundedCornerData& d);
void CalcOpenSection();
void AddElement(double b , double a);
void SetT(double t);
void TransformToNewFormat();
void Print();
void Underline(string a, char c='-', bool above=false);
void Underline(string a);
int NCells() const {return ns-nk+1;}
bool IsClosed() const {return NCells() >= 1;}
bool IsScheibenzug() const {return nk3 == 0;}
};
