lua-users home
lua-l archive

[Date Prev][Date Next][Thread Prev][Thread Next] [Date Index] [Thread Index]


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;}
};