Commit 4e55cde9 authored by thomas.forbriger's avatar thomas.forbriger Committed by thomas.forbriger
Browse files

proceeding

This is a legacy commit from before 2015-03-01.
It may be incomplete as well as inconsistent.
See COPYING.legacy and README.history for details.


SVN Path:     http://gpitrsvn.gpi.uni-karlsruhe.de/repos/TFSoftware/trunk
SVN Revision: 1252
SVN UUID:     67feda4a-a26e-11df-9d6e-31afc202ad0c
parent 2b2f41af
......@@ -3,7 +3,7 @@
*
* ----------------------------------------------------------------------------
*
* $Id: polymodel.cc,v 1.3 2002-12-30 18:08:29 forbrig Exp $
* $Id: polymodel.cc,v 1.4 2002-12-31 17:20:52 forbrig Exp $
* \author Thomas Forbriger
* \date 30/12/2002
*
......@@ -19,11 +19,12 @@
#define TF_POLYMODEL_CC_VERSION \
"TF_POLYMODEL_CC V1.0"
#define TF_POLYMODEL_CC_CVSID \
"$Id: polymodel.cc,v 1.3 2002-12-30 18:08:29 forbrig Exp $"
"$Id: polymodel.cc,v 1.4 2002-12-31 17:20:52 forbrig Exp $"
#include <gremlin1/polymodel.h>
#include <tfxx/error.h>
#include <aff/shaper.h>
#include <aff/slice.h>
#include <stdio.h>
#include <string>
#include <fstream>
......@@ -34,10 +35,16 @@ namespace gremlin1 {
//! define identifier
const char PolynomialModelFile::modversion2[];
//! define static
int PolynomialModelFile::norder;
//! define static
int PolynomialModelFile::nparameters;
//! default constructor
PolynomialModelFile::PolynomialModelFile():
Mpara(3,1,5), Mdepth(aff::Shaper(0,1)), Mfollow(1,5), Mnpol(1,5)
Mpara(norder,1,nparameters), Mdepth(aff::Shaper(0,1)),
Mfollow(1,nparameters), Mnpol(1,nparameters),
Mmeans(), Mmeans_are_valid(false)
{
Mpara=0.;
Mdepth=0.;
......@@ -63,17 +70,17 @@ namespace gremlin1 {
const char message[]="PolynomialModelFile: illegal index range";
try {
// check Mpara
TFXX_assert(((Mpara.f(0)==1)&&(Mpara.l(0)==3)), message);
TFXX_assert(((Mpara.f(0)==1)&&(Mpara.l(0)==norder)), message);
TFXX_assert(((Mpara.f(1)==1)&&(Mpara.l(1)>=Mpara.f(1))), message);
TFXX_assert(((Mpara.f(2)==1)&&(Mpara.l(2)==5)), message);
TFXX_assert(((Mpara.f(2)==1)&&(Mpara.l(2)==nparameters)), message);
// check Mdepth
TFXX_assert(((Mdepth.f(0)==0)&&(Mdepth.l(0)>Mdepth.f(0))), message);
// check Mfollow
TFXX_assert(((Mfollow.f(0)==1)&&(Mfollow.l(0)>=Mfollow.f(0))), message);
TFXX_assert(((Mfollow.f(1)==1)&&(Mfollow.l(1)==5)), message);
TFXX_assert(((Mfollow.f(1)==1)&&(Mfollow.l(1)==nparameters)), message);
// check Mnpol
TFXX_assert(((Mnpol.f(0)==1)&&(Mnpol.l(0)>=Mnpol.f(0))), message);
TFXX_assert(((Mnpol.f(1)==1)&&(Mnpol.l(1)==5)), message);
TFXX_assert(((Mnpol.f(1)==1)&&(Mnpol.l(1)==nparameters)), message);
// check for consistent number of sections
TFXX_assert(((Mpara.l(1)==Mdepth.l(0))&&(Mfollow.l(0)==Mdepth.l(0))&&
(Mnpol.l(0)==Mdepth.l(0))), message);
......@@ -91,6 +98,16 @@ namespace gremlin1 {
TFXX_assert((!(ishalfspace(i))), "illegal section index!");
return(Mdepth(valid_section_index(i)));
}
/*----------------------------------------------------------------------*/
//! return bottom depth of section
Tvalue PolynomialModelFile::top(const int& i) const
{
TFXX_assert((!(ishalfspace(i))), "illegal section index!");
return(Mdepth(valid_section_index(i)-1));
}
/*----------------------------------------------------------------------*/
//! return checked section index
......@@ -105,7 +122,7 @@ namespace gremlin1 {
//! return checked parameter index
int PolynomialModelFile::valid_parameter_index(const int& i) const
{
TFXX_assert(((i>0)&&(i<6)), "illegal parameter index!");
TFXX_assert(((i>0)&&(i<=nparameters)), "illegal parameter index!");
return(i);
}
......@@ -114,7 +131,7 @@ namespace gremlin1 {
//! return checked polynomial index
int PolynomialModelFile::valid_polynomial_index(const int& i) const
{
TFXX_assert(((i>0)&&(i<4)), "illegal polynomial index!");
TFXX_assert(((i>0)&&(i<=norder)), "illegal polynomial index!");
return(i);
}
......@@ -149,24 +166,6 @@ namespace gremlin1 {
/*----------------------------------------------------------------------*/
//! calculate parameter value at given depth
Tvalue PolynomialModelFile::value(const int& ipar,
const Tvalue& depth,
const int& isec) const
{
int is=valid_section_index(isec);
int ip=valid_parameter_index(ipar);
Tvalue retval=Mpara(1,is,ip);
if (!ishalfspace(is))
{
Tvalue x=polyarg(depth,is);
retval=Mpara(1,is,ip)+x*(Mpara(2,is,ip)+x*Mpara(3,is,ip));
}
return(retval);
}
/*----------------------------------------------------------------------*/
//! set unused polynomial coefficients to zero and check values
void PolynomialModelFile::clean_unused_and_check()
{
......@@ -228,6 +227,7 @@ namespace gremlin1 {
}
}
}
Mmeans_are_valid=false;
}
/*----------------------------------------------------------------------*/
......@@ -256,7 +256,7 @@ namespace gremlin1 {
os << " " << "<-- bottom of section";
os.width(5);
os << isec << " / polynomial expansion:" << std::endl;
for (int ipar=1; ipar<6; ++ipar)
for (int ipar=1; ipar<=nparameters; ++ipar)
{
os.width(12);
int npol=Mnpol(isec,ipar);
......@@ -264,9 +264,9 @@ namespace gremlin1 {
os << npol << " ";
}
os << std::endl << std::endl;
for (int ipol=1; ipol<4; ++ipol)
for (int ipol=1; ipol<=norder; ++ipol)
{
for (int ipar=1; ipar<6; ++ipar)
for (int ipar=1; ipar<=nparameters; ++ipar)
{
if (ipar<4)
{
......@@ -301,10 +301,10 @@ namespace gremlin1 {
is.ignore(maxskip,'\n');
// prepare arrays
try {
Mpara=Tarray(3,nsec+1,5);
Mpara=Tarray(norder,nsec+1,nparameters);
Mdepth=Tarray(aff::Shaper(0,nsec+1));
Mfollow=Tbarray(nsec+1,5);
Mnpol=Tiarray(nsec+1,5);
Mfollow=Tbarray(nsec+1,nparameters);
Mnpol=Tiarray(nsec+1,nparameters);
}
catch ( ... ) {
std::cerr << "PolynomialModelFile: could not create arrays for "
......@@ -324,7 +324,7 @@ namespace gremlin1 {
is.ignore(maxskip,'\n');
is >> Mdepth(isec);
is.ignore(maxskip,'\n');
for (int ipar=1; ipar<6; ++ipar)
for (int ipar=1; ipar<=nparameters; ++ipar)
{
is >> Mnpol(isec, ipar);
if (Mnpol(isec,ipar)<0)
......@@ -333,9 +333,9 @@ namespace gremlin1 {
Mfollow(isec,ipar)=true;
}
}
for (int ipol=1; ipol<4; ++ipol)
for (int ipol=1; ipol<=norder; ++ipol)
{
for (int ipar=1; ipar<6; ++ipar)
for (int ipar=1; ipar<=nparameters; ++ipar)
{
is >> Mpara(ipol, isec, ipar);
}
......@@ -343,7 +343,66 @@ namespace gremlin1 {
}
// finalize
clean_unused_and_check();
Mmeans_are_valid=false;
}
/*----------------------------------------------------------------------*/
//! calculate mean values for all sections
void PolynomialModelFile::calculate_means() const
{
if (!Mmeans_are_valid)
{
Mmeans=Tarray(norder,nsections(),nparameters);
for (int isec=Mmeans.f(1); isec<=Mmenas.l(1); ++isec)
{
for (int ipar=Mmeans.f(2); isec<=Mmeans.l(2); ++ipar)
{
Mmeans(3,isec,ipar)=2.*Mpara(3,isec,ipar);
Mmeans(2,isec,ipar)=Mpara(2,isec,ipar);
Mmeans(1,isec,ipar)=Mpara(1,isec,ipar)+
Mpara(3,isec,ipar)*thickness(isec)*thickness(isec)/3.;
}
}
Mmeans_are_valid=true;
}
}
/*----------------------------------------------------------------------*/
//! calculate mean values for depth range
/*
TCarray PolynomialModelFile::means(const Tvalue& top
const Tvalue& bottom) const
{
}
*/
/*----------------------------------------------------------------------*/
//! calculate derivative values at given depth
Tarray PolynomialModelFile::values(const Tvalue& depth,
const int& isec) const
{
int is=valid_section_index(isec);
Tarray result(norder,nparameters);
Tarray secpar(aff::slice(Mpara)()(is));
result.copyin(secpar);
if (!ishalfspace(is))
{
Tvalue x=polyarg(depth,is);
for (int ipar=result.f(1); ipar<=result.l(1); ++ipar)
{
result(3,ipar)=2.*secpar(3,ipar);
result(2,ipar)=secpar(2,ipar)+2.*secpar(3,ipar)*x;
result(1,ipar)=secpar(1,ipar)+x*(secpar(2,ipar)+x*secpar(3,ipar));
}
}
// NOTICE: returns reference (but result should have no reference to data
// members)
return(result);
}
} // namespace gremlin1
/* ----- END OF polymodel.cc ----- */
......@@ -3,7 +3,7 @@
*
* ----------------------------------------------------------------------------
*
* $Id: polymodel.h,v 1.4 2002-12-31 11:57:12 forbrig Exp $
* $Id: polymodel.h,v 1.5 2002-12-31 17:20:52 forbrig Exp $
* \author Thomas Forbriger
* \date 30/12/2002
*
......@@ -24,7 +24,7 @@
#define TF_POLYMODEL_H_VERSION \
"TF_POLYMODEL_H V1.0"
#define TF_POLYMODEL_H_CVSID \
"$Id: polymodel.h,v 1.4 2002-12-31 11:57:12 forbrig Exp $"
"$Id: polymodel.h,v 1.5 2002-12-31 17:20:52 forbrig Exp $"
#include<iostream>
#include<aff/array.h>
......@@ -53,6 +53,10 @@ namespace gremlin1 {
//! model file identifier
static const char modversion2[]="ModVersion 2";
//! number of polynomial orders
static const int norder=3;
//! number of physical parameters
static const int nparameters=5;
//! physical model parameters
enum Epara {
......@@ -82,6 +86,11 @@ namespace gremlin1 {
int nsections() const { return(Mpara.size(1)-1); }
//! lower interface of section \p i
Tvalue bottom(const int& i) const;
//! upper interface of section \p i
Tvalue top(const int& i) const;
//! thickness of section \p i
Tvalue thickness(const int& i) const
{ return(bottom(i)-top(i)); }
//! is this section index the half-space
bool ishalfspace(const int& i) const
{ return(i==Mdepth.l(0)); }
......@@ -103,16 +112,19 @@ namespace gremlin1 {
Tvalue value(const int& ipar,
const Tvalue& depth,
const int& isec) const;
//! return values for derivatives too at given depth
Tarray values(const Tvalue& depth) const
{ return(values(depth, secindex(depth))); }
//! return values for derivatives at given depth calculated from
//! polynomial for given section \p isec
Tarray values(const Tvalue& depth,
const int& isec) const;
//! return agrument of polynomial for given depth and section
Tvalue polyarg(const Tvalue& depth, const int& isec) const;
//! read access to parameters
TCarray parameters() const
{ return (aff::subarray(Mpara)()(nsections())()); }
//! read access to depth values
TCarray depth() const
{ return (aff::subarray(Mdepth)(1,nsections())); }
//! section means
TCarray sectionmeans() const
{ calculate_means(); return(Mmeans); }
private:
//! check for valid section index
int valid_section_index(const int& i) const;
......@@ -156,6 +168,17 @@ namespace gremlin1 {
* This is a matrix of size N x 5, where N-1 is the number of sections.
*/
Tiarray Mnpol;
/*----------------------------------------------------------------------*/
// some mutable data
// prepare polynomial means
void calculate_means() const;
// polynmial means within section
mutable Tarray Mmeans;
// are means valid
mutable bool Mmeans_are_valid;
}; // class PolynomialModelFile
/*----------------------------------------------------------------------*/
......
......@@ -3,7 +3,7 @@
*
* ----------------------------------------------------------------------------
*
* $Id: mocox.cc,v 1.2 2002-12-31 11:57:12 forbrig Exp $
* $Id: mocox.cc,v 1.3 2002-12-31 17:20:52 forbrig Exp $
* \author Thomas Forbriger
* \date 30/12/2002
*
......@@ -20,9 +20,10 @@
#define MOCOX_VERSION \
"MOCOX V1.1 convert gremlin1 model (modversion2)"
#define MOCOX_CVSID \
"$Id: mocox.cc,v 1.2 2002-12-31 11:57:12 forbrig Exp $"
"$Id: mocox.cc,v 1.3 2002-12-31 17:20:52 forbrig Exp $"
#include <iostream>
#include <cmath>
#include <fstream>
#include <gremlin1/polymodel.h>
#include <tfxx/commandline.h>
......@@ -48,8 +49,7 @@ typedef gremlin1::PolynomialModelFile::TCarray TCarray;
Tarray defnodes(const PolynomialModelFile& modelfile,
const bool& verbose)
{
if (verbose) cout << " define nodes:" << endl
<< " find appropriate number of nodes per section" << endl;
if (verbose) cout << " define nodes:" << endl;
// create local alias
int nsections=modelfile.nsections();
......@@ -79,19 +79,29 @@ Tarray defnodes(const PolynomialModelFile& modelfile,
*/
double scalefactor=0.01;
if (verbose) cout << " find appropriate number of nodes per section"
<< endl;
// array to store number of node layers per section
Tiarray nlay(nsections);
nlay=1;
int totalnlay=0;
// find number of layers per section
double bottom_of_upper=0.;
for (int i=1,i<=nsections; ++i)
{
if (verbose) cout << " sec #" << i <<":";
double d=depth(i)-bottom_of_upper;
if (verbose) cout << " d=" << d << "m";
// all second order coefficients in this section
TCarray c(aff::slice(para)(3)(i));
for (int j=c.f(0); j<=c.l(0); ++)
{
int n=1+int(d*sqrt(2.*abs(c(j))/(3.*scalefactor)));
nlay(i)= nlay(i)>n ? nlay(i) : n;
}
if (verbose) cout << " nlay=" << nlay(i);
totalnlay += nlay(i);
bottom_of_upper=depth(i);
}
} // defnodes
......
Markdown is supported
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