Subversion Repositories sifel

Compare Revisions

Ignore whitespace Rev 771 → Rev 772

/trunk/SIFEL/TRFEL/SRC/linbart3d.cpp
0,0 → 1,2618
/*
File: linbart.cpp
Author: Jaroslav Kruis, 31.3.2003
Purpose: onedimensional element with linear approximation functions
*/
 
#include "globalt.h"
#include "linbart.h"
#include "genfile.h"
#include "globmatt.h"
#include "loadcaset.h"
#include "loadelt.h"
#include "elemswitcht.h"
 
linbart3d::linbart3d (void)
{
long i;
// number of nodes on element
nne=2;
// geometric problem dimension (3D)
ncomp=3;
// number of end nodes
nen=2;
// the number of edges
ned=1;
// the number of nodes on one edge
nned=0;
 
// number of transported variables
ntm=Tp->ntm;
 
 
dofe = new long* [ntm];
nip = new long* [ntm];
intordkm = new long* [ntm];
intordcm = new long* [ntm];
ordering = new long* [ntm];
for (i=0;i<ntm;i++){
dofe[i] = new long [ntm];
nip[i] = new long [ntm];
intordkm[i] = new long [ntm];
intordcm[i] = new long [ntm];
ordering[i] = new long [nne];
}
switch (Tp->tmatt){
case onemedium:{
ordering[0][0]=1; ordering[0][1]=2;
dofe[0][0]=2; intordkm[0][0]=3; intordcm[0][0]=3; nip[0][0]=6;
ndofe=2; napfun=1;
break;
}
case twomediacoup:{
ordering[0][0]=1; ordering[0][1]=3;
ordering[1][0]=2; ordering[1][1]=4;
 
intordkm[0][0]=1; intordkm[0][1]=1; intordkm[1][0]=1; intordkm[1][1]=1;
intordcm[0][0]=2; intordcm[0][1]=2; intordcm[1][0]=2; intordcm[1][1]=2;
if (Tp->savemode==0){
nip[0][0]=3; nip[0][1]=3; nip[1][0]=3; nip[1][1]=3;
}
if (Tp->savemode==1){
nip[0][0]=3; nip[0][1]=0; nip[1][0]=0; nip[1][1]=0;
}
dofe[0][0]=2; dofe[0][1]=2; dofe[1][0]=2; dofe[1][1]=2;
ndofe=4; napfun=2;
break;
}
case threemediacoup:{
ordering[0][0]=1; ordering[0][1]=4;
ordering[1][0]=2; ordering[1][1]=5;
ordering[2][0]=3; ordering[2][1]=6;
 
intordkm[0][0]=1; intordkm[0][1]=1; intordkm[0][2]=1;
intordkm[1][0]=1; intordkm[1][1]=1; intordkm[1][2]=1;
intordkm[2][0]=1; intordkm[2][1]=1; intordkm[2][2]=1;
 
intordcm[0][0]=2; intordcm[0][1]=2; intordcm[0][2]=2;
intordcm[1][0]=2; intordcm[1][1]=2; intordcm[1][2]=2;
intordcm[2][0]=2; intordcm[2][1]=2; intordcm[2][2]=2;
 
if (Tp->savemode==0){
nip[0][0]=3; nip[0][1]=3; nip[0][2]=3;
nip[1][0]=3; nip[1][1]=3; nip[1][2]=3;
nip[2][0]=3; nip[2][1]=3; nip[2][2]=3;
}
if (Tp->savemode==1){
nip[0][0]=3; nip[0][1]=0; nip[0][2]=0;
nip[1][0]=0; nip[1][1]=0; nip[1][2]=0;
nip[2][0]=0; nip[2][1]=0; nip[2][2]=0;
}
dofe[0][0]=2; dofe[0][1]=2; dofe[0][2]=2;
dofe[1][0]=2; dofe[1][1]=2; dofe[1][2]=2;
dofe[2][0]=2; dofe[2][1]=2; dofe[2][2]=2;
 
ndofe=6; napfun=3;
break;
}
case fourmediacoup:{
ordering[0][0]=1; ordering[0][1]=5;
ordering[1][0]=2; ordering[1][1]=6;
ordering[2][0]=3; ordering[2][1]=7;
ordering[3][0]=4; ordering[3][1]=8;
 
intordkm[0][0]=1; intordkm[0][1]=1; intordkm[0][2]=1; intordkm[0][3]=1;
intordkm[1][0]=1; intordkm[1][1]=1; intordkm[1][2]=1; intordkm[1][3]=1;
intordkm[2][0]=1; intordkm[2][1]=1; intordkm[2][2]=1; intordkm[2][3]=1;
intordkm[3][0]=1; intordkm[3][1]=1; intordkm[3][2]=1; intordkm[3][3]=1;
 
intordcm[0][0]=2; intordcm[0][1]=2; intordcm[0][2]=2; intordcm[0][3]=2;
intordcm[1][0]=2; intordcm[1][1]=2; intordcm[1][2]=2; intordcm[1][3]=2;
intordcm[2][0]=2; intordcm[2][1]=2; intordcm[2][2]=2; intordcm[2][3]=2;
intordcm[3][0]=2; intordcm[3][1]=2; intordcm[3][2]=2; intordcm[3][3]=2;
 
if (Tp->savemode==0){
nip[0][0]=3; nip[0][1]=3; nip[0][2]=3; nip[0][3]=3;
nip[1][0]=3; nip[1][1]=3; nip[1][2]=3; nip[1][3]=3;
nip[2][0]=3; nip[2][1]=3; nip[2][2]=3; nip[2][3]=3;
nip[3][0]=3; nip[3][1]=3; nip[3][2]=3; nip[3][3]=3;
}
if (Tp->savemode==1){
nip[0][0]=3; nip[0][1]=0; nip[0][2]=0; nip[0][3]=0;
nip[1][0]=0; nip[1][1]=0; nip[1][2]=0; nip[1][3]=0;
nip[2][0]=0; nip[2][1]=0; nip[2][2]=0; nip[2][3]=0;
nip[3][0]=0; nip[3][1]=0; nip[3][2]=0; nip[3][3]=0;
}
dofe[0][0]=2; dofe[0][1]=2; dofe[0][2]=2; dofe[0][3]=2;
dofe[1][0]=2; dofe[1][1]=2; dofe[1][2]=2; dofe[1][3]=2;
dofe[2][0]=2; dofe[2][1]=2; dofe[2][2]=2; dofe[2][3]=2;
dofe[3][0]=2; dofe[3][1]=2; dofe[3][2]=2; dofe[3][3]=2;
 
ndofe=8; napfun=4;
break;
}
default:{
print_err("unknown type of transported matter is required",__FILE__,__LINE__,__func__);
}
}
}
 
linbart3d::~linbart3d (void)
{
long i;
 
for (i=0;i<ntm;i++){
delete [] dofe[i];
delete [] nip[i];
delete [] intordkm[i];
delete [] intordcm[i];
delete [] ordering[i];
}
delete [] dofe;
delete [] nip;
delete [] intordkm;
delete [] intordcm;
delete [] ordering;
}
 
 
 
/**
Function computes direction %vector.
@param s - direction %vector
@param x,y,z - vectors of nodal coordinates
@retval Function returns direction %vector s.
 
TKr 06/08/2021 accroding to JK
*/
void linbart3d::dirvect (vector &s,vector &x,vector &y,vector &z)
{
double l;
 
s[0]=x[1]-x[0];
s[1]=y[1]-y[0];
s[2]=z[1]-z[0];
l=sqrt(s[0]*s[0]+s[1]*s[1]+s[2]*s[2]);
if (l<Tp->zero){
print_err("zero norm of direction vector",__FILE__,__LINE__,__func__);
}
s[0]/=l;
s[1]/=l;
s[2]/=l;
}
 
 
/**
Function returns coordinates of nodes at local element coordinate system.
 
@param x - %vector of nodal x-coordinates in global coordinate system
@param y - %vector of nodal y-coordinates in global coordinate system
@param z - %vector of nodal z-coordinates in global coordinate system
@param lx - %vector of nodal coordinates at local element coordinate system
 
@retval Function returns %vector of nodal coordinates at parameter lx.
 
TKr 06/08/2021 accroding to JK
*/
void linbart3d::giveloccoord(vector &x, vector &y, vector &z, vector &lx)
{
double l;
 
lx[0] = 0.0;
l = sqr(x[1] - x[0]) + sqr(y[1] - y[0]) + sqr(z[1] - z[0]);
l = sqrt(l);
lx[1] = l;
}
 
 
/**
function assembles code numbers for one medium
@param cn - code numbers
@param ri - row index
JK
*/
void linbart3d::codnum (long *cn,long ri)
{
long i;
for (i=0;i<nne;i++){
cn[i]=ordering[ri][i];
}
}
 
/**
function approximates function defined by nodal values
 
@param xi - natural coordinate on element
@param nodval - %vector of nodal values
TKr 06/08/2021 accroding to JK
*/
double linbart3d::approx (double xi,vector &nodval)
{
double f;
vector bf(ASTCKVEC(nne));
bf_lin_1d (bf.a,xi);
scprd (bf,nodval,f);
return f;
}
 
/**
function computes values in integration points from nodal values
 
@param eid - element id
 
TKr 06/08/2021 accroding to JK
*/
void linbart3d::intpointval (long eid)
{
long i,k,ii,jj,ipp;
double xi,val;
vector r(ASTCKVEC(ndofe)),t(ASTCKVEC(nne)),gp,w;
elemvalues(eid, r);
for (k=0;k<ntm;k++){
// nodal values of a single variable
for (i=0;i<dofe[k][k];i++){
t[i]=r[ordering[k][i]-1];
}
for (ii=0;ii<Tp->ntm;ii++){
for (jj=0;jj<Tp->ntm;jj++){
// integration points for the conductivity matrix
reallocv (RSTCKVEC(intordkm[ii][jj],gp));
reallocv (RSTCKVEC(intordkm[ii][jj],w));
gauss_points (gp.a,w.a,intordkm[ii][jj]);
ipp=Tt->elements[eid].ipp[ii][jj];
for (i=0;i<intordkm[ii][jj];i++){
xi=gp[i];
val = approx (xi,t);
Tm->ip[ipp].av[k]=val;
ipp++;
}
// integration points for the capacity matrix
reallocv (RSTCKVEC(intordcm[ii][jj],gp));
reallocv (RSTCKVEC(intordcm[ii][jj],w));
gauss_points (gp.a,w.a,intordcm[ii][jj]);
for (i=0;i<intordcm[ii][jj];i++){
xi=gp[i];
val = approx (xi,t);
Tm->ip[ipp].av[k]=val;
ipp++;
}
if (Tp->savemode==1) break;
}
if (Tp->savemode==1) break;
}
}
}
 
 
/**
function computes values in integration points from nodal values
this function is used for arbitrary variable, not only for
variables used as unknowns in the problem
@param eid - element id
@param nodval - nodal values
@param ipval - value at integration points
@retval Function returns approximated values at integration points in the %vector ipval.
 
TKr 06/08/2021 accroding to JK
*/
void linbart3d::intpointval (long eid,vector &nodval,vector &ipval)
{
long i,kk,ii,jj,ipp;
double xi;
vector gp,w;
kk = 0;
for (ii=0;ii<Tp->ntm;ii++){
for (jj=0;jj<Tp->ntm;jj++){
// interpolation to integration points of conductivity matrix
reallocv (RSTCKVEC(intordkm[ii][jj],gp));
reallocv (RSTCKVEC(intordkm[ii][jj],w));
gauss_points (gp.a,w.a,intordkm[ii][jj]);
ipp=Tt->elements[eid].ipp[ii][jj];
for (i=0;i<intordkm[ii][jj];i++){
xi=gp[i];
// value in integration point
ipval[kk] = approx (xi,nodval);
ipp++;
kk++;
}
// interpolation to integration points of capacity matrix
reallocv (RSTCKVEC(intordcm[ii][jj],gp));
reallocv (RSTCKVEC(intordcm[ii][jj],w));
gauss_points (gp.a,w.a,intordcm[ii][jj]);
for (i=0;i<intordcm[ii][jj];i++){
xi=gp[i];
// value in integration point
ipval[kk] = approx (xi,nodval);
 
ipp++;
kk++;
}
if (Tp->savemode==1) break;
}
if (Tp->savemode==1) break;
}
}
 
 
/**
function computes initial values in integration points from initial nodal values
 
@param eid - element id
 
TKr 06/08/2021 accroding to TKo
*/
void linbart3d::initintpointval (long eid)
{
long i,k,ii,jj,ipp,ndofn,cndofn;
double xi,val;
ivector enod(ASTCKIVEC(nne));
vector rr, r(ASTCKVEC(ndofe)), t(ASTCKVEC(nne)), gp, w;
Tt->give_elemnodes (eid,enod);
 
for(i=0, cndofn=0; i<nne; i++)
{
ndofn = Tt->give_ndofn(enod[i]);
// make reference rr to the vector of nodal values on element
makerefv(rr, r.a+cndofn, ndofn);
// get initial nodal values and store them in the vector r with the help of reference rr
initnodval2(enod[i], rr);
cndofn += ndofn;
}
for (k=0;k<ntm;k++){
// nodal values of a single variable
for (i=0;i<dofe[k][k];i++){
t[i]=r[ordering[k][i]-1];
}
for (ii=0;ii<Tp->ntm;ii++){
for (jj=0;jj<Tp->ntm;jj++){
// integration points for the conductivity matrix
reallocv (RSTCKVEC(intordkm[ii][jj],gp));
reallocv (RSTCKVEC(intordkm[ii][jj],w));
gauss_points (gp.a,w.a,intordkm[ii][jj]);
ipp=Tt->elements[eid].ipp[ii][jj];
for (i=0;i<intordkm[ii][jj];i++){
xi=gp[i];
val = approx (xi,t);
Tm->ip[ipp].av[k]=val;
ipp++;
}
// integration points for the capacity matrix
reallocv (RSTCKVEC(intordcm[ii][jj],gp));
reallocv (RSTCKVEC(intordcm[ii][jj],w));
gauss_points (gp.a,w.a,intordcm[ii][jj]);
for (i=0;i<intordcm[ii][jj];i++){
xi=gp[i];
val = approx (xi,t);
Tm->ip[ipp].av[k]=val;
ipp++;
}
if (Tp->savemode==1) break;
}
if (Tp->savemode==1) break;
}
}
}
 
 
/**
function computes values in integration points from nodal values
 
@param eid - element id
 
TKr 06/08/2021 accroding to JK
*/
void linbart3d::intpointgrad (long eid)
{
long i,ii,jj,k,ipp;
double xi,jac;
vector x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),r(ASTCKVEC(ndofe)),t(ASTCKVEC(nne));
vector gp,w,grad(ASTCKVEC(ncomp));
matrix gm(ncomp,nne);
Tt->give_node_coord3d (x,y,z,eid);
elemvalues(eid, r);
 
for (k=0;k<Tp->ntm;k++){
// nodal values of a single variable
for (i=0;i<dofe[k][k];i++){
t[i]=r[ordering[k][i]-1];
}
for (ii=0;ii<Tp->ntm;ii++){
for (jj=0;jj<Tp->ntm;jj++){
// integration points for the conductivity matrix
reallocv (RSTCKVEC(intordkm[ii][jj],gp));
reallocv (RSTCKVEC(intordkm[ii][jj],w));
gauss_points (gp.a,w.a,intordkm[ii][jj]);
ipp=Tt->elements[eid].ipp[ii][jj];
for (i=0;i<intordkm[ii][jj];i++){
xi=gp[i];
// matrix of gradients
grad_matrix (gm,x,y,z,xi,jac);
mxv (gm,t,grad);
Tm->storegrad (k,ipp,grad);
ipp++;
}
// integration points for the capacity matrix
reallocv (RSTCKVEC(intordcm[ii][jj],gp));
reallocv (RSTCKVEC(intordcm[ii][jj],w));
gauss_points (gp.a,w.a,intordcm[ii][jj]);
for (i=0;i<intordcm[ii][jj];i++){
xi=gp[i];
// matrix of gradients
grad_matrix (gm,x,y,z,xi,jac);
mxv (gm,t,grad);
Tm->storegrad (k,ipp,grad);
ipp++;
}
if (Tp->savemode==1) break;
}
if (Tp->savemode==1) break;
}
}
}
 
/**
function approximates nodal values of array other to integration points
@param eid - element id
TKr 06/08/2021 accroding to JK
*/
void linbart3d::intpointother (long eid)
{
long i, k, ii, jj, ipp, ncompo, nodid;
double xi, val;
ivector nodes(nne);
vector t(ASTCKVEC(nne)), r, gp, w;
// nodes of required element
Tt->give_elemnodes (eid,nodes);
// first node number
nodid=nodes[0];
 
// number of components
ncompo=Tt->nodes[nodid].ncompother;
reallocv(RSTCKVEC(ncompo*nne, r));
// nodal values of array other
nodalotherval (nodes, r);
for (k=0;k<ncompo;k++){
// nodal values of a single variable
for (i=0;i<nne;i++){
t[i]=r[i*ncompo+k];
}
for (ii=0;ii<Tp->ntm;ii++){
for (jj=0;jj<Tp->ntm;jj++){
// integration points for the conductivity matrix
 
reallocv (intordkm[ii][jj],gp);
reallocv (intordkm[ii][jj],w);
gauss_points (gp.a,w.a,intordkm[ii][jj]);
ipp=Tt->elements[eid].ipp[ii][jj];
for (i=0;i<intordkm[ii][jj];i++){
xi=gp[i];
val = approx (xi,t);
Tm->ip[ipp].other[k]=val;
ipp++;
}
// integration points for the capacity matrix
 
reallocv (intordcm[ii][jj],gp);
reallocv (intordcm[ii][jj],w);
gauss_points (gp.a,w.a,intordcm[ii][jj]);
for (i=0;i<intordcm[ii][jj];i++){
xi=gp[i];
val = approx (xi,t);
Tm->ip[ipp].other[k]=val;
ipp++;
}
if (Tp->savemode==1) break;
}
if (Tp->savemode==1) break;
}
}
}
 
/**
function assembles approximation or test functions into a %matrix
@param w - %matrix of test functions
@param xi - natural coordinate
@param v - velocity of advection
TKr 06/08/2021 accroding to JK
*/
void linbart3d::btf_matrix (matrix &n,double xi,double v)
{
if (Tp->advect==0)
bf_matrix (n,xi);
if (Tp->advect==1)
tf_matrix (n,xi,v);
}
 
 
/**
function assembles %matrix of base functions
@param n - %matrix of base functions
@param xi - natural coordinate
TKr 06/08/2021 accroding to JK
*/
void linbart3d::bf_matrix (matrix &n,double xi)
{
nullm (n);
bf_lin_1d (n.a,xi);
}
 
/**
function assembles %matrix of test functions
@param w - %matrix of test functions
@param xi - natural coordinate
@param v - velocity of advection
TKr 06/08/2021 accroding to JK
*/
void linbart3d::tf_matrix (matrix &w,double xi,double v)
{
double alpha=1.0;
matrix n(1,dofe[0][0]);
nullm (n);
bf_lin_1d (n.a,xi);
 
double h=1.2500000000e-03;
double nor=1.507500e-07;
double pe=nor*h/2.0/9.52e-12;
alpha=(exp(pe)+exp(0.0-pe))/(exp(pe)-exp(0.0-pe))-1.0/pe;
//alpha=1.0;
 
if (v<0.0){
w[0][0]=n[0][0]+alpha/2.0;
w[0][1]=n[0][1]-alpha/2.0;
}else{
w[0][0]=n[0][0]-alpha/2.0;
w[0][1]=n[0][1]+alpha/2.0;
}
}
 
 
/**
function assembles gradient of %matrix of base functions
@param gm - gradient %matrix
@param x - array containing node coordinates
@param xi - natural coordinate
@param jac - jacobian
TKr 06/08/2021 accroding to JK
*/
void linbart3d::grad_matrix (matrix &gm,vector &x,vector &y,vector &z,double xi,double &jac)
{
long i;
vector s(ASTCKVEC(3)),lx(ASTCKVEC(nne)),dx(ASTCKVEC(nne));
 
// direction vector
dirvect (s,x,y,z);
// local coordinates of nodes
giveloccoord (x,y,z,lx);
// derivatives of the approximation functions
derivatives_1d (dx,jac,lx,xi);
for (i=0;i<nne;i++){
gm[0][i]=dx[i]*s[0];
gm[1][i]=dx[i]*s[1];
gm[2][i]=dx[i]*s[2];
}
}
 
/**
function computes conductivity %matrix of 1D problems for one transported matter
finite element with linear approximation functions
@param lcid - number of load case
@param eid - number of element
@param ri,ci - row and column indices of the computed block in the resulting %matrix
@param km - conductivity %matrix
TKr 06/08/2021 accroding to JK
*/
void linbart3d::conductivity_matrix (long lcid,long eid,long ri,long ci,matrix &km)
{
long i,ii;
double area,xi,ww,jac;
ivector nodes(nne);
vector x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),w(intordkm[ri][ci]),gp(intordkm[ri][ci]),a(nne);
matrix gm(ncomp,dofe[ri][ci]),d;
matrix n(1,dofe[ri][ci]);
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// node coordinates
Tt->give_node_coord3d (x,y,z,eid);
// coordinates and weights of integration points
gauss_points (gp.a,w.a,intordkm[ri][ci]);
nullm (km);
if (Tp->savemode==0)
ii=Tt->elements[eid].ipp[ri][ci];
if (Tp->savemode==1)
ii=Tt->elements[eid].ipp[0][0];
for (i=0;i<intordkm[ri][ci];i++){
xi=gp[i]; ww=w[i];
// matrix of gradients
grad_matrix (gm,x,y,z,xi,jac);
// matrix of conductivity of the material
reallocm(ncomp,ncomp,d);
Tm->matcond (d,ii,ri,ci);
 
// area of cross section
area = approx (xi,a);
jac*=area*ww;
// contribution to the conductivity matrix of the element
bdbj (km.a,gm.a,d.a,jac,gm.m,gm.n);
//convective terms
reallocm (1,ncomp,d);
Tm->matcond2 (d,ii,ri,ci);
btf_matrix (n, xi,d[0][0]);
bdbjac (km, n, d, gm, jac);
ii++;
}
if ((Tt->elements[eid].transi[lcid]==3) || (Tt->elements[eid].transi[lcid]==4)){
// additional matrix due to transmission
transmission_matrix (lcid,eid,ri,ci,km);
}
}
 
/**
function computes capacity %matrix of 1D problems for one transported matter
finite element with bilinear approximation functions
@param eid - number of element
@param ri,ci - row and column indices of the computed block in the resulting %matrix
@param cm - capacity %matrix
TKr 06/08/2021 accroding to JK
*/
void linbart3d::capacity_matrix (long eid,long ri,long ci,matrix &cm)
{
long i,ii;
double jac,xi,ww,area,rho,c;
ivector nodes(nne);
vector lx(ASTCKVEC(nne)),s(ASTCKVEC(3));
vector x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),w(intordcm[ri][ci]),gp(intordcm[ri][ci]),a(nne);
matrix n(1,dofe[ri][ci]);
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
Tc->give_densitye (eid,rho);
// node coordinates
Tt->give_node_coord3d (x,y,z,eid);
 
// local coordinates of nodes
giveloccoord(x,y,z,lx);
 
// coordinates and weights of integration points
gauss_points (gp.a,w.a,intordcm[ri][ci]);
nullm (cm);
 
if (Tp->savemode==0)
ii=Tt->elements[eid].ipp[ri][ci]+intordkm[ri][ci];
if (Tp->savemode==1)
ii=Tt->elements[eid].ipp[0][0]+intordkm[0][0];
 
for (i=0;i<intordcm[ri][ci];i++){
xi=gp[i]; ww=w[i];
bf_matrix (n,xi);
// area of cross section
area = approx (xi,a);
c=Tm->capcoeff (ii,ri,ci);
jac_1d (jac,lx,xi);
jac*=area*ww*rho*c;
nnj (cm.a,n.a,jac,n.m,n.n);
ii++;
}
}
 
 
/**
function computes source %vector of one matter on one element
\int_{Omega} N^T N d Omega . s
@param sv - source %vector of one matter
@param nodval - array of nodal values
@param eid - element id
@param ri,ci - row and column indices of the block (ri must be equal to ci)
TKr 06/08/2021 accroding to JK
*/
void linbart3d::quantity_source_vector (vector &sv,vector &nodval,long eid,long ri,long ci)
{
long i;
double jac,xi,ww,area;
ivector nodes(nne);
vector lx(ASTCKVEC(nne)),x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),w(intordcm[ri][ci]),gp(intordcm[ri][ci]),v(dofe[ri][ci]),a(nne);
matrix n(1,dofe[ri][ci]),nm(dofe[ri][ci],dofe[ri][ci]);
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// node coordinates
Tt->give_node_coord3d (x,y,z,eid);
 
// local coordinates of nodes
giveloccoord(x,y,z,lx);
 
// coordinates and weights of integration points
gauss_points (gp.a,w.a,intordcm[ri][ci]);
nullm (nm);
for (i=0;i<intordcm[ri][ci];i++){
xi=gp[i]; ww=w[i];
bf_matrix (n,xi);
jac_1d (jac,lx,xi);
// area of cross section
area = approx (xi,a);
 
jac*=ww*area;
nnj (nm.a,n.a,jac,n.m,n.n);
}
mxv (nm,nodval,v); addv (sv,v,sv);
}
 
 
 
/**
function computes transmission complement to the conductivity %matrix for one matter
\int_{Gamma_3} N^T c_{tr} N dGamma
@param lcid - number of load case
@param eid - number of element
@param ri,ci - row and column indices
@param km - part of the conductivity %matrix
TKr 06/08/2021 accroding to JK
*/
void linbart3d::transmission_matrix (long lcid,long eid,long ri,long ci,matrix &km)
{
long i,neleml,ii;
bocontypet bc[2]; // array for boundary condition indicators
 
double xi,area;
double new_trc;//added
long nn;//added
 
ivector nodes(nne);
vector trc(ASTCKVEC(2)),a(ASTCKVEC(nne));
matrix n(1,dofe[ri][ci]);
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
for (i=0;i<Tb->lc[lcid].neb;i++){
if (Tb->lc[lcid].elemload[i].eid==eid){
neleml=i;
break;
}
}
 
// indicators of boundary conditions
Tb->lc[lcid].elemload[neleml].give_bc (bc);
// transmission coefficients
Tb->lc[lcid].elemload[neleml].give_trc (Tp->time,Tb->lc[lcid].nodval,lcid,ci,trc);
if (Tp->savemode==0)
ii=Tt->elements[eid].ipp[ri][ci];
if (Tp->savemode==1)
ii=Tt->elements[eid].ipp[0][0];
 
if (bc[0]==5 || bc[0]>10){
xi=-1.0;
bf_matrix (n,xi);
/***********************************************///added
nn = nodes[0];
Tm->transmission_transcoeff(new_trc,trc[0],ri,ci,nn,bc[0],ii);
//printf("new_trc = %e\n",new_trc);
 
trc[0]=new_trc;
/***********************************************///added
// area of cross section
area = approx (xi,a);
nnj (km.a,n.a,trc[0]*area,n.m,n.n);
}
 
if (bc[1]==5 || bc[1]>10){
xi=1.0;
bf_matrix (n,xi);
 
/***********************************************///added
nn = nodes[1];
Tm->transmission_transcoeff(new_trc,trc[1],ri,ci,nn,bc[1],ii);
//printf("new_trc = %e\n",new_trc);
trc[1]=new_trc;
/***********************************************///added
// area of cross section
area = approx (xi,a);
 
nnj (km.a,n.a,trc[1]*area,n.m,n.n);
}
}
 
 
 
/**
function computes contributions to the transmission %vector
\int_{Gamma_3} N^T c_{tr} N dGamma * nodal_external_value
@param tmv - transmission %vector of one matter
@param lcid - load case id (corresponds to the row index)
@param eid - element id
@param leid - loaded element id
@param cid - component id (corresponds to the column index)
TKr 06/08/2021 accroding to JK
*/
void linbart3d::transmission_vector (vector &tmv,long lcid,long eid,long leid,long cid)
{
long ii,nid;
bocontypet bc[2]; // array for boundary condition indicators
double xi,area,new_nodval,tr;
ivector nodes(nne);
vector trc(ASTCKVEC(2)),trcn(ASTCKVEC(nne)),nodval(ASTCKVEC(2)),av(ASTCKVEC(dofe[lcid][cid])),v(ASTCKVEC(dofe[lcid][cid])),a(nne);
vector trr(ASTCKVEC(2));
matrix n(1,dofe[lcid][cid]),km(dofe[lcid][cid],dofe[lcid][cid]);
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// indicators of boundary conditions
Tb->lc[lcid].elemload[leid].give_bc (bc);
// nodal values of exterior variable
Tb->lc[lcid].elemload[leid].give_external_nodval (lcid,cid,nodval);
// transmission coefficients
Tb->lc[lcid].elemload[leid].give_trc (Tp->time,Tb->lc[lcid].nodval,lcid,cid,trc);
// transmission/radiation coefficients
Tb->lc[lcid].elemload[leid].give_trr (Tp->time,Tb->lc[lcid].nodval,trr);
 
// integration point number
if (Tp->savemode==0)
ii=Tt->elements[eid].ipp[lcid][cid];
if (Tp->savemode==1)
ii=Tt->elements[eid].ipp[0][0];
 
if (bc[0]==5 || bc[0]>10){
nullm (km);
xi=-1.0;
bf_matrix (n,xi);
 
// node id
nid = nodes[0];
tr = 0.0;
if(bc[0] == 90)
tr = trr[0]/trc[0];
Tm->transmission_nodval (new_nodval,nodval[0],tr,lcid,cid,nid,bc[0],ii);
av[0]=new_nodval; av[1]=0.0;
// area of cross section
area = approx (xi,a);
//fprintf (stdout,"\n linbart3d trc[0] %le",trc[0]);
//fprintf (stdout,"\n linbart3d val[0] %le",av[0]);
 
nnj (km.a,n.a,trc[0]*area,n.m,n.n);
mxv (km,av,v); addv (tmv,v,tmv);
}
if (bc[1]==5 || bc[1]>10){
nullm (km);
xi=1.0;
bf_matrix (n,xi);
 
/***********************************************///added
nid = nodes[1];
tr = 0.0;
if(bc[0] == 90)
tr = trr[1]/trc[1];
Tm->transmission_nodval(new_nodval,nodval[1],tr,lcid,cid,nid,bc[1],ii);
av[0]=0.0; av[1]=new_nodval;
/***********************************************/
// area of cross section
area = approx (xi,a);
//fprintf (stdout,"\n linbart3d trc[1] %le",trc[1]);
//fprintf (stdout,"\n linbart3d val[1] %le",av[1]);
nnj (km.a,n.a,trc[1]*area,n.m,n.n);
 
mxv (km,av,v); addv (tmv,v,tmv);
}
}
 
 
 
/**
function computes contribution to the convection %vector
 
\int_{Gamma_2} N^T N dGamma * nodal_flux_values
@param f - convection %vector of one matter
@param lcid - load case id
@param eid - element id
@param leid - loaded element id
@param ri,ci - row and column indices of the block (ri must be equal to ci)
 
TKr 06/08/2021 accroding to JK
*/
void linbart3d::convection_vector (vector &f,long lcid,long eid,long leid,long ri,long ci)
{
bocontypet bc[2]; // array for boundary condition indicators
double xi,area;
ivector nodes(nne);
vector nodval(ASTCKVEC(2)),av(ASTCKVEC(dofe[ri][ci])),v(ASTCKVEC(dofe[ri][ci])),a(ASTCKVEC(nne));
matrix n(1,dofe[ri][ci]),km(dofe[ri][ci],dofe[ri][ci]);
 
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// indicators of boundary conditions
Tb->lc[lcid].elemload[leid].give_bc (bc);
// nodal values of exterior variable
Tb->lc[lcid].elemload[leid].give_nodval (lcid,nodval);
 
 
if (bc[0]==2 || bc[0]==3 || bc[0]==5){
nullm (km);
xi=-1.0;
 
bf_matrix (n,xi);
// area of cross section
area = approx (xi,a);
nnj (km.a,n.a,area,n.m,n.n);
 
av[0]=nodval[0]; av[1]=0.0;
mxv (km,av,v); addv (f,v,f);
}
if (bc[1]==2 || bc[1]==3 || bc[1]==5){
nullm (km);
xi=1.0;
bf_matrix (n,xi);
// area of cross section
area = approx (xi,a);
nnj (km.a,n.a,area,n.m,n.n);
 
av[0]=0.0; av[1]=nodval[1];
mxv (km,av,v); addv (f,v,f);
}
}
 
/**
function computes internal fluxes of 1D problems for one transported matter
finite element with linear approximation functions
@param lcid - number of load case
@param eid - number of element
@param ifl - %vector of internal fluxes
TKr 06/08/2021 accroding to JK
*/
void linbart3d::internal_fluxes (long lcid,long eid,vector &ifl)
{
long i,ipp;
double jac,area,xi;
ivector nodes(nne);
vector x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),w(ASTCKVEC(intordkm[lcid][lcid])),gp(ASTCKVEC(intordkm[lcid][lcid])),fl(ASTCKVEC(ncomp)),contr(ASTCKVEC(dofe[lcid][lcid])),a(ASTCKVEC(nne));
matrix gm(ncomp,dofe[lcid][lcid]),d(ncomp,ncomp);
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// node coordinates
Tt->give_node_coord3d (x,y,z,eid);
// coordinates and weights of integration points
gauss_points (gp.a,w.a,intordkm[lcid][lcid]);
if (Tp->savemode==0)
ipp=Tt->elements[eid].ipp[lcid][lcid];
if (Tp->savemode==1)
ipp=Tt->elements[eid].ipp[0][0];
 
 
for (i=0;i<intordkm[lcid][lcid];i++){
xi=gp[i];
Tm->computenlfluxes (lcid,ipp);
Tm->givefluxes (lcid,ipp,fl);
 
grad_matrix (gm,x,y,z,gp[i],jac);
mtxv (gm,fl,contr);
// area of cross section
area = approx (xi,a);
cmulv (area*jac*w[i],contr);
addv (contr,ifl,ifl);
ipp++;
}
}
 
/**
function computes energy accumulated in an element
finite element with linear approximation functions
@param lcid - number of load case
@param eid - number of element
@param ifl - %vector of internal fluxes
TKr 06/08/2021 accroding to JK
*/
void linbart3d::accumul_energy (long lcid,long eid,vector &acum)
{
long i,ipp;
double jac,area,xi,ae,v;
ivector nodes(nne);
vector lx(ASTCKVEC(nne)),x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),w(ASTCKVEC(intordcm[lcid][lcid])),gp(ASTCKVEC(intordcm[lcid][lcid])),contr(ASTCKVEC(dofe[lcid][lcid])),a(ASTCKVEC(nne));
matrix n(1,dofe[lcid][lcid]);
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// node coordinates
Tt->give_node_coord3d (x,y,z,eid);
 
// local coordinates of nodes
giveloccoord(x,y,z,lx);
 
// coordinates and weights of integration points
gauss_points (gp.a,w.a,intordcm[lcid][lcid]);
if (Tp->savemode==0)
ipp=Tt->elements[eid].ipp[lcid][lcid]+intordkm[lcid][lcid];
if (Tp->savemode==1)
ipp=Tt->elements[eid].ipp[0][0]+intordkm[0][0];
 
 
for (i=0;i<intordcm[lcid][lcid];i++){
xi=gp[i];
Tm->computenlcapacities (lcid,ipp);
// velocity has to be obtained from somewhere
v=0.0;
 
btf_matrix (n,xi,v);
contr[0]=n[0][0];
contr[1]=n[0][1];
 
ae = Tm->ip[ipp].cap[lcid];
// area of cross section
area = approx (xi,a);
jac_1d (jac,lx,xi);
cmulv (ae*area*jac*w[i],contr);
addv (contr,acum,acum);
ipp++;
}
}
 
 
/**
function computes advection %matrix of 1D problems for one transported matter
finite element with linear approximation functions
@param lcid - number of load case
@param eid - number of element
@param ri,ci - row and column indices of the computed block in the resulting %matrix
@param hm - advection %matrix
TKr 06/08/2021 accroding to JK
*/
void linbart3d::advection_matrix (long lcid,long eid,long ri,long ci,matrix &hm)
{
long i,ii;
double area,xi,ww,jac;
ivector nodes(nne);
vector x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),w(ASTCKVEC(intordkm[ri][ci])),gp(ASTCKVEC(intordkm[ri][ci])),a(ASTCKVEC(nne));
matrix v(1,ncomp);
matrix gm(ncomp,dofe[ri][ci]), n(1,dofe[ri][ci]);
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// node coordinates
Tt->give_node_coord3d (x,y,z,eid);
// coordinates and weights of integration points
gauss_points (gp.a,w.a,intordkm[ri][ci]);
nullm (hm);
if (Tp->savemode==0)
ii=Tt->elements[eid].ipp[ri][ci];
if (Tp->savemode==1)
ii=Tt->elements[eid].ipp[0][0];
for (i=0;i<intordkm[ri][ci];i++){
xi=gp[i]; ww=w[i];
// matrix of approximation functions
bf_matrix (n,xi);
// matrix of gradients
grad_matrix (gm,x,y,z,xi,jac);
// vector of velocity
Tm->matcond2(v,ii,ri,ci);
// area of cross section
area = approx (xi,a);
jac*=area*ww;
// contribution to the advection matrix of the element
bdbjac (hm,n,v,gm,jac);
ii++;
}
}
 
 
 
/**
function assembles resulting element conductivity %matrix
 
@param eid - element id
@param lcid - load case id
@param km - resulting conductivity %matrix of one element
 
TKr 06/08/2021 accroding to JK
*/
void linbart3d::res_conductivity_matrix (long eid,long lcid,matrix &km)
{
long i,j,*rcn,*ccn;
matrix lkm;
for (i=0;i<ntm;i++){
for (j=0;j<ntm;j++){
rcn = new long [dofe[i][j]];
ccn = new long [dofe[i][j]];
reallocm (dofe[i][j],dofe[i][j],lkm);
conductivity_matrix (i,eid,i,j,lkm);
codnum (rcn,i);
codnum (ccn,j);
mat_localize (km,lkm,rcn,ccn);
delete [] rcn; delete [] ccn;
}
}
}
 
 
 
/**
function computes contributions to the right-hand side - volume integral
\int_{Omega} B^T D dOmega
@param tmv - transmission %vector of one matter
@param lcid - load case id
@param eid - element id
@param ri,ci - row and column indices of the block (ri must be equal to ci)
TKr 06/08/2021
*/
void linbart3d::volume_rhs_vector (long lcid,long eid,long ri,long ci,vector &vrhs)
{
long i,ii;
double area,xi,ww,jac;
ivector nodes(nne);
vector x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),w(ASTCKVEC(intordkm[ri][ci])),gp(ASTCKVEC(intordkm[ri][ci])),a(ASTCKVEC(nne));
matrix gm(ncomp,dofe[ri][ci]),d;
matrix km(dofe[ri][ci],1);
 
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// node coordinates
Tt->give_node_coord3d (x,y,z,eid);
// coordinates and weights of integration points
gauss_points (gp.a,w.a,intordkm[ri][ci]);
nullm (km);
if (Tp->savemode==0)
ii=Tt->elements[eid].ipp[ri][ci];
if (Tp->savemode==1)
ii=Tt->elements[eid].ipp[0][0];
for (i=0;i<intordkm[ri][ci];i++){
xi=gp[i]; ww=w[i];
// matrix of gradients
grad_matrix (gm,x,y,z,xi,jac);
// matrix of conductivity of the material
reallocm(ncomp,1,d);
 
Tm->volume_rhs (d,ii,ri,ci,ncomp);
// area of cross section
area = approx (xi,a);
 
jac*=area*ww;
// contribution to the volume_rhs integral of the element
nnjac (km,gm,d,jac);
ii++;
}
 
for (i=0;i<vrhs.n;i++){
vrhs[i] = km[i][0];
}
 
}
 
 
 
/**
function computes contributions to the right-hand side - volume integral of the second type
\int_{Omega} N^T const dOmega
@param vrhs - volume right-hand side %vector of one matter
@param lcid - load case id
@param eid - element id
@param ri,ci - row and column indices of the block (ri must be equal to ci)
TKr 06/08/2021
*/
void linbart3d::volume_rhs_vector2 (long lcid,long eid,long ri,long ci,vector &vrhs)
{
long i,ii,ipp;
double jac,xi,w1,area,c;
ivector nodes(nne);
vector lx(ASTCKVEC(nne)),x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),w(ASTCKVEC(intordcm[ri][ci])),gp(ASTCKVEC(intordcm[ri][ci])),a(ASTCKVEC(nne));
matrix n(1,dofe[ri][ci]);
vector nn(dofe[ri][ci]);
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// node coordinates
Tt->give_node_coord3d (x,y,z,eid);
 
// local coordinates of nodes
giveloccoord(x,y,z,lx);
 
// coordinates and weights of integration points
gauss_points (gp.a,w.a,intordcm[ri][ci]);
nullv (vrhs);
if (Tp->savemode==0)
ipp=Tt->elements[eid].ipp[ri][ci];
if (Tp->savemode==1)
ipp=Tt->elements[eid].ipp[0][0];
for (i=0;i<intordcm[ri][ci];i++){
xi=gp[i]; w1=w[i];
jac_1d (jac,lx,xi);
bf_matrix (n,xi);
c=Tm->volume_rhs2 (ipp,ri,ci);
area = approx (xi,a);
jac*=w1*area*c;
for (ii=0;ii<n.n;ii++){
nn[ii] = n[0][ii];
}
addmultv(vrhs,nn,jac);
ipp++;
}
}
 
 
 
 
/**
function assembles resulting element volume right-hand side
 
@param eid - element id
@param lcid - load case id
@param f - resulting volume right-hand side %vector of one element
 
TKr 06/08/2021
*/
void linbart3d::res_volume_rhs_vector (vector &f,long eid,long lcid)
{
long i,*cn;
vector lf;
 
for (i=0;i<ntm;i++){
cn = new long [dofe[i][i]];
reallocv (dofe[i][i],lf);
codnum (cn,i);
volume_rhs_vector (i,eid,i,i,lf);
locglob (f.a,lf.a,cn,dofe[i][i]);
delete [] cn;
}
}
 
 
/**
function assembles resulting element volume right-hand side of the second type
 
@param eid - element id
@param lcid - load case id
@param f - resulting volume right-hand side %vector of one element
 
TKr 06/08/2021
*/
void linbart3d::res_volume_rhs_vector2 (vector &f,long eid,long lcid)
{
long i,*cn;
vector lf;
 
for (i=0;i<ntm;i++){
cn = new long [dofe[i][i]];
reallocv (dofe[i][i],lf);
codnum (cn,i);
volume_rhs_vector2 (i,eid,i,i,lf);
locglob (f.a,lf.a,cn,dofe[i][i]);
delete [] cn;
}
}
 
 
 
 
/**
function assembles resulting element capacity %matrix
@param eid - element id
@param cm - resulting capacity %matrix of one element
 
TKr 06/08/2021 accroding to JK
*/
void linbart3d::res_capacity_matrix (long eid,matrix &cm)
{
long i,j,*rcn,*ccn;
matrix lcm;
for (i=0;i<ntm;i++){
for (j=0;j<ntm;j++){
rcn = new long [dofe[i][j]];
ccn = new long [dofe[i][j]];
reallocm (dofe[i][j],dofe[i][j],lcm);
capacity_matrix (eid,i,j,lcm);
 
// diagonalization of capacity matrix on one element
if (Tp->diagcap == 1){
diagonalization (lcm);
}
codnum (rcn,i);
codnum (ccn,j);
mat_localize (cm,lcm,rcn,ccn);
delete [] rcn; delete [] ccn;
}
}
}
 
/**
function assembles resulting element convection %vector
 
@param f - resulting convection %vector of one element
@param lcid - load case id
@param eid - element id
@param leid - element id
TKr 06/08/2021 accroding to JK
*/
void linbart3d::res_convection_vector (vector &f,long lcid,long eid,long leid)
{
long *cn;
vector lf;
 
// transi[lcid]==2 - element contains boundary with prescribed flux
// transi[lcid]==4 - element contains boundary with prescribed transmission and prescribed flux
if ((Tt->elements[eid].transi[lcid]==2)||(Tt->elements[eid].transi[lcid]==4)){
// array for code numbers
cn = new long [dofe[lcid][lcid]];
// array for nodal values of one matter
reallocv (dofe[lcid][lcid],lf);
 
convection_vector (lf,lcid,eid,leid,lcid,lcid);
 
// code numbers
codnum (cn,lcid);
// localization to element vector
locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
delete [] cn;
cmulv (-1.0,f,f);
}
}
 
/**
function assembles resulting element transmission %vector
 
@param f - resulting transmission %vector of one element
@param lcid - load case id
@param eid - element id
@param leid - loaded element id
TKr 06/08/2021 accroding to JK
*/
void linbart3d::res_transmission_vector (vector &f,long lcid,long eid,long leid)
{
long *cn,i;
vector lf;
cn = new long [dofe[lcid][lcid]];
// code numbers
codnum (cn,lcid);
reallocv (dofe[lcid][lcid],lf);
// transi[i]==3 - element contains boundary with prescribed transmission
// transi[i]==4 - element contains boundary with prescribed transmission and prescribed flux
if ((Tt->elements[eid].transi[lcid]==3)||(Tt->elements[eid].transi[lcid]==4)){
for (i=0;i<ntm;i++){
nullv (lf.a,lf.n);
 
transmission_vector (lf,lcid,eid,leid,i);
locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
}
}
delete [] cn;
}
 
/**
function assembles resulting element source %vector
 
@param sv - resulting source %vector of one element
@param lcid - load case id
@param eid - element id
TKr 06/08/2021 accroding to JK
*/
void linbart3d::res_quantity_source_vector (vector &sv,vector &nodval,long lcid,long eid)
{
long *cn;
vector lsv;
 
cn = new long [dofe[lcid][lcid]];
reallocv (dofe[lcid][lcid],lsv);
codnum (cn,lcid);
quantity_source_vector (lsv,nodval,eid,lcid,lcid);
locglob (sv.a,lsv.a,cn,dofe[lcid][lcid]);
delete [] cn;
}
 
 
 
/**
function assembles resulting element internal fluxes vector
 
@param eid - element id
@param elemif - resulting internal fluxes %vector of one element
TKr 06/08/2021 accroding to JK
*/
void linbart3d::res_internal_fluxes (long eid,vector &elemif)
{
long i;
vector lif, tdnv, nodval;
ivector cn;
matrix cm;
 
for (i=0; i<ntm; i++){
reallocv(RSTCKIVEC(dofe[i][i], cn));
reallocv(RSTCKVEC(dofe[i][i], lif));
codnum(cn.a, i);
internal_fluxes(i, eid, lif);
locglob(elemif.a, lif.a, cn.a, dofe[i][i]);
}
 
reallocv(RSTCKVEC(ndofe, tdnv));
reallocv(RSTCKVEC(ndofe, lif));
reallocm(RSTCKMAT(ndofe, ndofe, cm));
nodalderivatives(eid, tdnv);
res_capacity_matrix(eid, cm);
mxv(cm, tdnv, lif);
subv(lif, elemif, elemif);
if (Tt->elements[eid].source==1){
reallocv(RSTCKVEC(nne, nodval));
nullv(nodval);
nullv(lif);
Tb->lc[0].sourcenodalvalues(eid, nodval);
source_vector(0, eid, nodval, lif);
addv(elemif, lif, elemif);
}
}
 
 
/**
function computes nodal energies for one transported matter
finite element with linear approximation functions
@param lcid - number of load case
@param eid - number of element
@param nener - %vector of nodal energy
TKr 06/08/2021 accroding to JK
*/
void linbart3d::nodal_energies (long lcid,long eid,vector &nener)
{
long i,ipp;
double jac,area,xi;
ivector nodes(nne);
vector x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),w(ASTCKVEC(intordkm[lcid][lcid])),gp(ASTCKVEC(intordkm[lcid][lcid])),fl(ASTCKVEC(ncomp)),contr(ASTCKVEC(dofe[lcid][lcid])),a(ASTCKVEC(nne));
matrix gm(ncomp,dofe[lcid][lcid]),d(ncomp,ncomp);
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// node coordinates
Tt->give_node_coord3d (x,y,z,eid);
// coordinates and weights of integration points
gauss_points (gp.a,w.a,intordkm[lcid][lcid]);
if (Tp->savemode==0)
ipp=Tt->elements[eid].ipp[lcid][lcid];
if (Tp->savemode==1)
ipp=Tt->elements[eid].ipp[0][0];
 
 
for (i=0;i<intordkm[lcid][lcid];i++){
xi=gp[i];
Tm->computenlfluxes (lcid,ipp);
Tm->givefluxes (lcid,ipp,fl);
 
grad_matrix (gm,x,y,z,gp[i],jac);
mtxv (gm,fl,contr);
// area of cross section
area = approx (xi,a);
cmulv (area*jac*w[i],contr);
addv (contr,nener,nener);
ipp++;
}
}
 
 
/**
function assembles resulting nodal energies in nodes on element
 
@param eid - element id
@param elemne - resulting %vector of nodal energies
TKr 06/08/2021 accroding to JK
*/
void linbart3d::res_nodal_energy (long eid,vector &elemne,double dt)
{
long i;
vector nener, tdnv, nodval;
ivector cn;
matrix cm;
 
for (i=0; i<ntm; i++){
reallocv(RSTCKIVEC(dofe[i][i], cn));
reallocv(RSTCKVEC(dofe[i][i], nener));
codnum(cn.a, i);
internal_fluxes(i, eid, nener);
locglob(elemne.a, nener.a, cn.a, dofe[i][i]);
}
cmulv (dt,elemne);
reallocv(RSTCKVEC(ndofe, tdnv));
reallocv(RSTCKVEC(ndofe, nener));
reallocm(RSTCKMAT(ndofe, ndofe, cm));
//nodalderivatives(eid, tdnv);
//res_capacity_matrix(eid, cm);
//mxv(cm, tdnv, nener);
 
for (i=0; i<ntm; i++){
reallocv(RSTCKIVEC(dofe[i][i], cn));
reallocv(RSTCKVEC(dofe[i][i], nener));
codnum(cn.a, i);
accumul_energy (i,eid,nener);
cmulv(-1.0,nener);
locglob(elemne.a, nener.a, cn.a, dofe[i][i]);
}
if (Tt->elements[eid].source==1){
reallocv(RSTCKVEC(nne, nodval));
nullv(nodval);
nullv(nener);
Tb->lc[0].sourcenodalvalues(eid, nodval);
source_vector(0, eid, nodval, nener);
addv(elemne, nener, elemne);
}
}
 
 
 
/**
function assembles resulting element advection %matrix
 
@param eid - element id
@param lcid - load case id
@param hm - resulting advection %matrix of one element
 
TKr 06/08/2021 accroding to JK
*/
void linbart3d::res_advection_matrix (long eid,long lcid,matrix &hm)
{
long i,j,*rcn,*ccn;
matrix lhm;
for (i=0;i<ntm;i++){
for (j=0;j<ntm;j++){
rcn = new long [dofe[i][j]];
ccn = new long [dofe[i][j]];
reallocm (dofe[i][j],dofe[i][j],lhm);
advection_matrix (i,eid,i,j,lhm);
codnum (rcn,i);
codnum (ccn,j);
mat_localize (hm,lhm,rcn,ccn);
delete [] rcn; delete [] ccn;
}
}
}
 
 
/**
function computes element quantity integral
the quantity is stored in nodes
 
@param eid - element id
@param nodval - %vector of quantity nodal values
 
@retval f - element quantity integral
 
TKr 06/08/2021
*/
double linbart3d::total_integral (long eid,vector &nodval)
{
long i;
double area,xi,ww,jac,value,f;
ivector nodes(nne);
vector lx(ASTCKVEC(nne)),x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),w(ASTCKVEC(2)),gp(ASTCKVEC(2)),a(ASTCKVEC(nne));
 
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// node coordinates
Tt->give_node_coord3d (x,y,z,eid);
 
// local coordinates of nodes
giveloccoord(x,y,z,lx);
gauss_points (gp.a,w.a,2);
f = 0.0;
 
for (i=0;i<2;i++){
xi=gp[i]; ww=w[i];
jac_1d (jac,lx,xi);
 
value = approx (xi,nodval);
// area of cross section
area = approx (xi,a);
 
jac*=area*ww*value;
f = f + jac;
}
return(f);
}
 
/**
function computes element quantity integral
the quantity is stored in integration points
 
@param eid - element id
@param varid - id of variable required
@retval f - element quantity integral
 
TKr 06/08/2021 accroding to JK
*/
double linbart3d::total_integral_ip (long eid,long varid)
{
long i,ii;
double area,xi,ww,jac,value,f;
ivector nodes(nne);
vector lx(ASTCKVEC(nne)),x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne)),w(ASTCKVEC(2)),gp(ASTCKVEC(2)),a(ASTCKVEC(nne));
 
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// node coordinates
Tt->give_node_coord3d (x,y,z,eid);
 
// local coordinates of nodes
giveloccoord(x,y,z,lx);
gauss_points (gp.a,w.a,2);
f = 0.0;
long ri=0;
long ci=0;
if (Tp->savemode==0)
ii=Tt->elements[eid].ipp[ri][ci];
if (Tp->savemode==1)
ii=Tt->elements[eid].ipp[0][0];
 
for (i=0;i<2;i++){
xi=gp[i]; ww=w[i];
jac_1d (jac,lx,xi);
 
value = Tm->ip[ii].eqother[varid];
// area of cross section
area = approx (xi,a);
 
jac*=area*ww*value;
f = f + jac;
}
return f;
}
 
/**************************************************************///added
/**
function computes contributions to the boundary flux from transmission %vector
\int_{Gamma_3} N^T c_{tr} N dGamma * nodal_external_value
 
@param tmv - %vector of boundary fluxes
@param lcid - load case id
@param eid - element id
@param leid - element id
@param ri,ci - row and column indices of the computed block in the resulting %matrix
TKr 06/08/2021
*/
void linbart3d::boundary_flux (vector &tmv,long lcid,long eid,long leid,long ri,long ci)
{
long ii;
double xi,area;
double new_nodval,tr;
long nn;
bocontypet bc[2];
ivector nodes(nne);
vector trc(ASTCKVEC(2)),trcn(ASTCKVEC(nne)),nodval(ASTCKVEC(2)),av(ASTCKVEC(dofe[ri][ci])),v(ASTCKVEC(dofe[ri][ci])),a(ASTCKVEC(nne));
vector trr(ASTCKVEC(2));
matrix n(1,dofe[ri][ci]),km(dofe[ri][ci],dofe[ri][ci]);
// element nodes
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
// indicators of boundary conditions
Tb->lc[lcid].elemload[leid].give_bc (bc);
// nodal values of exterior variable
Tb->lc[lcid].elemload[leid].give_nodval (lcid,nodval);
// transmission coefficients
Tb->lc[lcid].elemload[leid].give_trc (Tp->time,Tb->lc[lcid].nodval,lcid,ci,trc);
// transmission/radiation coefficients
Tb->lc[lcid].elemload[leid].give_trr (Tp->time,Tb->lc[lcid].nodval,trr);
 
if (Tp->savemode==0)
ii=Tt->elements[eid].ipp[ri][ci];
if (Tp->savemode==1)
ii=Tt->elements[eid].ipp[0][0];
 
if (bc[0]==5 || bc[0]>10){
nullm (km);
xi=-1.0;
bf_matrix (n,xi);
 
nn = nodes[0];
tr = 0.0;
if(bc[0] == 90)
tr = trr[0]/trc[0];
Tm->transmission_flux(new_nodval,nodval[0],tr,ri,ci,nn,bc[0],ii);
av[0]=new_nodval; av[1]=0.0;
// area of cross section
area = approx (xi,a);
nnj (km.a,n.a,trc[0]*area,n.m,n.n);
mxv (km,av,v); addv (tmv,v,tmv);
}
if (bc[1]==5 || bc[1]>10){
nullm (km);
xi=1.0;
bf_matrix (n,xi);
 
nn = nodes[1];
tr = 0.0;
if(bc[0] == 90)
tr = trr[1]/trc[1];
Tm->transmission_flux(new_nodval,nodval[1],tr,ri,ci,nn,bc[1],ii);
av[0]=0.0; av[1]=new_nodval;
// area of cross section
area = approx (xi,a);
nnj (km.a,n.a,trc[1]*area,n.m,n.n);
 
mxv (km,av,v); addv (tmv,v,tmv);
}
}
 
 
 
/**
function assembles resulting element boundary flux %vector
 
@param f - resulting boundary flux %vector of one element
@param lcid - load case id
@param eid - element id
@param leid - element id
TKr 06/08/2021
*/
void linbart3d::res_boundary_flux (vector &f,long lcid,long eid,long leid)
{
long *cn,i;
vector lf;
 
cn = new long [dofe[lcid][lcid]];
codnum (cn,lcid);
//coupled b.c.
for (i=0;i<ntm;i++){
reallocv (dofe[lcid][lcid],lf);
if ((Tt->elements[eid].transi[i]==3)||(Tt->elements[eid].transi[i]==4))
boundary_flux (lf,i,eid,leid,lcid,i);
locglob (f.a,lf.a,cn,dofe[lcid][lcid]);
}
delete [] cn;
}
 
 
/**
function computes correct fluxes at integration points on element
 
@param eid - element id
TKr 06/08/2021
*/
void linbart3d::intpointflux (long eid)
{
long i,ii,jj,k,ipp;
for (k=0;k<Tp->ntm;k++){
for (ii=0;ii<Tp->ntm;ii++){
for (jj=0;jj<Tp->ntm;jj++){
ipp=Tt->elements[eid].ipp[ii][jj];
for (i=0;i<intordkm[ii][jj];i++){
// computation of correct fluxes
if (Tp->fluxcomp==1)
Tm->computenlfluxes (k,ipp);
ipp++;
}
for (i=0;i<intordcm[ii][jj];i++){
// computation of correct fluxes
if (Tp->fluxcomp==1)
Tm->computenlfluxes (k,ipp);
ipp++;
}
if (Tp->savemode==1) break;
}
if (Tp->savemode==1) break;
}
}
}
 
/**
function computes gradients in nodes of element
 
@param eid - element id
TKr 06/08/2021
*/
void linbart3d::nod_grads_ip (long eid)
{
long i,j,k,ipp;
ivector ipnum(nne),nod(nne);
vector grad(ASTCKVEC(ncomp));
// numbers of integration points closest to nodes
// (function is from the file GEFEL/ordering.cpp)
ipp=Tt->elements[eid].ipp[0][0];
nodip_bar (ipp,intordkm[0][0],ipnum);
 
// node numbers of the element
Tt->give_elemnodes (eid,nod);
for (i=0;i<nne;i++){
for (k=0;k<Tp->ntm;k++){
// strains at the closest integration point
Tm->givegrad (k,ipnum[i],grad);
// storage of strains to the node
j=nod[i];
Tt->nodes[j].storegrad (k,grad);
}
}
}
 
 
/**
function computes fluxes in nodes of element
 
@param eid - element id
TKr 06/08/2021
*/
void linbart3d::nod_fluxes_ip (long eid)
{
long i,j,k,ipp;
ivector ipnum(nne),nod(nne);
vector flux(ASTCKVEC(ncomp));
// numbers of integration points closest to nodes
// (function is from the file GEFEL/ordering.cpp)
ipp=Tt->elements[eid].ipp[0][0];
nodip_bar (ipp,intordkm[0][0],ipnum);
// node numbers of the element
Tt->give_elemnodes (eid,nod);
intpointflux (eid);
 
for (i=0;i<nne;i++){
for (k=0;k<Tp->ntm;k++){
// strains at the closest integration point
Tm->givefluxes (k,ipnum[i],flux);
// storage of strains to the node
j=nod[i];
Tt->nodes[j].storeflux (k,flux);
}
}
}
 
 
 
/**
function computes other values directly in nodes of element
 
@param lcid - load case id
@param eid - element id
@param ri - row index
@param ci - column index
 
TKr 06/08/2021
*/
void linbart3d::nod_others_comp (long lcid,long eid,long ri,long ci)
{
long i,j,k,ipp,ncompother,ii;
vector h,r(ASTCKVEC(ndofe));
ivector nod(ASTCKIVEC(nne));
double other;
Tt->give_elemnodes (eid,nod);
elemvalues(eid, r);
ncompother = Tm->givencompother();
reallocv (RSTCKVEC(ntm,h));
 
ii = 0;
for (i=0;i<nne;i++){
if (Tp->savemode==0)
ipp=Tt->elements[eid].ipp[ri][ci];
if (Tp->savemode==1)
ipp=Tt->elements[eid].ipp[0][0];
for (j=0;j<ntm;j++){
h[j] = r[ii];
ii++;
}
for(k=0;k<ncompother;k++){
other = Tm->givecompother(k,ipp,h.a);
Tt->nodes[nod[i]].storeother (k,other);
}
}
}
 
 
 
/**
Function returns transport (non-mechanical) quantities at nodes of element.
The values of selected quantity is copied from the closest integration points
to element nodes.
 
@param eid - element id
@param nodval - %vector of nodal values
@param ntq - type of non-mechanical quantity
@return The function does not return anything.
TKr 06/08/2021 accroding to TKo
*/
void linbart3d::transq_nodval (long eid,vector &nodval,nonmechquant nmq)
{
long i,ipid;
ivector ipnum(nne);
// numbers of integration points closest to nodes
// (function is from the file GEFEL/ordering.cpp)
ipid=Tt->elements[eid].ipp[0][0];
nodip_bar (ipid,intordkm[0][0],ipnum);
 
for (i=0;i<nne;i++)
{
//copy transport (non-mechanical) quantity from closest int. point
nodval[i] = Tm->givetransq(nmq, ipnum[i]);
}
}
 
 
 
/**
Function computes transport (non-mechanical) quantities in nodes of element.
 
@param eid - element id
@param nodval - %vector of nodal values of all required quantities, i.e.,
nodal value of i-th quantity in j-th node is given by nodval[i*ncnv+j] where ncnv
is the number of calculated nodes on eid-th element.
@param ncne - number of computed nodes on element (only first ncne of nodes is calculated)
@param nq - number of required transport quantities
@param qt - array of types of required non-mechanical, i.e. transport quantities
@return The function does not return anything.
TKr 06/08/2021 accroding to TKo
*/
void linbart3d::transq_nodval_comp (long eid,vector &nodval, long ncne, long nq, nonmechquant *qt)
{
long i,j,ipid;
long ncompo;
ivector ipnum(nne), enod(nne);
intpointst ipb; // backup of the first integration point of element
vector ipav;
// numbers of integration points closest to nodes
// (function is from the file GEFEL/ordering.cpp)
ipid=Tt->elements[eid].ipp[0][0];
 
// element nodes
Tt->give_elemnodes(eid, enod);
 
// numbers of integration points closest to element nodes
// (function is from the file GEFEL/ordering.cpp)
nodip_bar (ipid,intordkm[0][0],ipnum);
 
// calculate nodal gradients - actually, only one rarely used material model depends on gradients
// nod_grads_comp(eid);
 
// store original content of the first integration point on element because
// it becomes working int. point for nodal values calculations on the given element
ipb.copy(Tm->ip[ipid], ntm, 1);
 
// The first integration point will be used for computation of nodal values temporarily
// then the original content of the first integration point will be restored
for (i=0;i<ncne;i++)
{
makerefv(ipav, Tm->ip[ipid].av, Tp->ntm);
// store nodal values to the first (working) integration point
nodalval (enod[i], ipav);
// store gradients to the first (working) integration point
// for(j=0; j<ntm; j++)
// Tm->ip[ipid].storegrad(j, Tt->nodes[enod[i]].gradient[j]);
ncompo = Tm->ip[ipnum[i]].ncompeqother;
if (ncompo){
// take eqother values for the given node from the closest integration point
Tm->storeeqother(ipid, 0, ncompo, Tm->ip[ipnum[i]].eqother);
}
ncompo = Tm->ip[ipnum[i]].ncompother;
if (ncompo){
// take eqother values for the given node from the closest integration point
Tm->storeother(ipid, 0, ncompo, Tm->ip[ipnum[i]].other);
}
Tm->mat_aux_values(ipid);
 
//give calculated transport quantity from the first int. point
for (j=0; j<nq; j++)
nodval[j*ncne+i] = Tm->givetransq(qt[j], ipid);
}
// restore original integration point content of other/eqother arrays
Tm->ip[ipid].copy(ipb, ntm, 0);
}
 
 
 
/**
Function returns initial values of transport (non-mechanical) quantities at nodes of element.
The values of selected quantity is copied from the closest integration points
to element nodes.
 
@param eid - element id
@param nodval - %vector of nodal values
@param ntq - type of non-mechanical quantity
@return The function does not return anything.
TKr 06/08/2021 accroding to TKo
*/
void linbart3d::transq_init_nodval (long eid,vector &nodval,nonmechquant nmq)
{
long i,ipid;
ivector ipnum(ASTCKIVEC(nne)), enod(ASTCKIVEC(nne));
vector ipav;
// numbers of integration points closest to nodes
// (function is from the file GEFEL/ordering.cpp)
ipid=Tt->elements[eid].ipp[0][0];
nodip_bar (ipid,intordkm[0][0],ipnum);
// element nodes
Tt->give_elemnodes(eid, enod);
 
for (i=0;i<nne;i++)
{
// create reference vector to the int. point av array
makerefv(ipav, Tm->ip[ipnum[i]].av, Tp->ntm);
// store initial nodal values to the integration point
initnodval2(enod[i], ipav);
//copy transport (non-mechanical) quantity from closest int. point
nodval[i] = Tm->givetransq(nmq, ipnum[i]);
}
}
 
 
 
/**
Function computes initial transport (non-mechanical) quantities in nodes of element.
 
@param eid - element id
@param nodval - %vector of nodal values of all required quantities, i.e.,
nodal value of i-th quantity in j-th node is given by nodval[i*ncnv+j] where ncnv
is the number of calculated nodes on eid-th element.
@param ncne - number of computed nodes on element (only first ncne of nodes is calculated)
@param nq - number of required transport quantities
@param qt - array of types of required non-mechanical, i.e. transport quantities
@return The function does not return anything.
TKr 06/08/2021 accroding to TKo
*/
void linbart3d::transq_init_nodval_comp (long eid,vector &nodval, long ncne, long nq, nonmechquant *qt)
{
long i,j,ipid;
long ncompo;
ivector ipnum(ASTCKIVEC(nne)), enod(ASTCKIVEC(nne));
intpointst ipb; // backup of the first integration point of element
vector ipav;
// numbers of integration points closest to nodes
// (function is from the file GEFEL/ordering.cpp)
ipid=Tt->elements[eid].ipp[0][0];
 
// element nodes
Tt->give_elemnodes(eid, enod);
 
// numbers of integration points closest to element nodes
// (function is from the file GEFEL/ordering.cpp)
nodip_bar (ipid,intordkm[0][0],ipnum);
 
// calculate nodal gradients - actually, only one rarely used material model depends on gradients
// nod_grads_comp(eid);
 
// store original content of the first integration point on element because
// it becomes working int. point for nodal values calculations on the given element
ipb.copy(Tm->ip[ipid], ntm, 1);
 
// The first integration point will be used for computation of nodal values temporarily
// then the original content of the first integration point will be restored
for (i=0;i<ncne;i++)
{
// create reference vector to the int. point av array
makerefv(ipav, Tm->ip[ipid].av, Tp->ntm);
// store initial nodal values to the first (working) integration point
initnodval2 (enod[i], ipav);
// store gradients to the first (working) integration point
// for(j=0; j<ntm; j++)
// Tm->ip[ipid].storegrad(j, Tt->nodes[enod[i]].gradient[j]);
ncompo = Tm->ip[ipnum[i]].ncompeqother;
if (ncompo){
// take eqother values for the given node from the closest integration point
Tm->storeeqother(ipid, 0, ncompo, Tm->ip[ipnum[i]].eqother);
}
ncompo = Tm->ip[ipnum[i]].ncompother;
if (ncompo){
// take eqother values for the given node from the closest integration point
Tm->storeother(ipid, 0, ncompo, Tm->ip[ipnum[i]].other);
}
Tm->mat_aux_values(ipid);
 
//give calculated transport quantity from the first int. point
for (j=0; j<nq; j++)
nodval[j*ncne+i] = Tm->givetransq(qt[j], ipid);
}
// restore original integration point content of other/eqother arrays
Tm->ip[ipid].copy(ipb, ntm, 0);
}
 
 
 
/**
Function assembles global coordinates of integration points.
@param eid - element id
@param ipp - integration point pointer
@param ri - row index of the block of integration points /according to transported media/ (input)
@param ci - column index of the block of integration points /according to transported media/ (input)
@param coord - array containing coordinates of integration points (output)
@retval Function returns coordinates of integration points in the %vector coord.
@retval 0
 
TKr 06/08/2021 accroding to TKo
*/
long linbart3d::ipcoord (long eid, long ipp, long ri, long ci, vector &coord)
{
long i, ii;
double xi;
vector x(ASTCKVEC(nne)), y(ASTCKVEC(nne)), w, gp;
 
Tt->give_node_coord2d(x, y, eid);
// integration points for the conductivity matrix
reallocv(RSTCKVEC(intordkm[ri][ci], gp));
reallocv(RSTCKVEC(intordkm[ri][ci], w));
gauss_points(gp.a, w.a, intordkm[ri][ci]);
ii = Tt->elements[eid].ipp[ri][ci];
for (i=0; i<intordkm[ri][ci]; i++)
{
if (ii == ipp)
{
xi=gp[i];
coord[0] = approx(xi, x);
coord[1] = approx(xi, y);
coord[2] = 0.0;
return 0;
}
ii++;
}
// integration points for the capacity matrix
reallocv(RSTCKVEC(intordcm[ri][ci], gp));
reallocv(RSTCKVEC(intordcm[ri][ci], w));
gauss_points (gp.a, w.a, intordcm[ri][ci]);
for (i=0; i<intordcm[ri][ci]; i++)
{
if (ii == ipp)
{
xi = gp[i];
coord[0] = approx(xi, x);
coord[1] = approx(xi, y);
coord[2] = 0.0;
return 0;
}
ii++;
}
 
return 1;
}
 
 
 
/**
Function computes natural coordinates of integration point.
@param eid - element id (input)
@param ipp - integration point pointer (input)
@param ncoord - array containing coordinates of integration points (output)
@retval coord - function returns coordinates of integration points in the %vector ncoord.
@retval 0 - on success
@retval 1 - integration point with ipp was not found
 
TKr 06/08/2021 accroding to TKo
*/
long linbart3d::ipncoord (long eid, long ipp, vector &ncoord)
{
long i, ri, ci, ii;
double xi;
vector w, gp;
 
for (ri=0;ri<ntm;ri++)
{
for (ci=0;ci<ntm;ci++)
{
// integration points for the conductivity matrix
reallocv(RSTCKVEC(intordkm[ri][ci], gp));
reallocv(RSTCKVEC(intordkm[ri][ci], w));
gauss_points(gp.a, w.a, intordkm[ri][ci]);
ii = Tt->elements[eid].ipp[ri][ci];
for (i=0; i<intordkm[ri][ci]; i++)
{
if (ii == ipp)
{
xi=gp[i];
ncoord[0] = xi;
ncoord[1] = 0.0;
ncoord[2] = 0.0;
return 0;
}
ii++;
}
// integration points for the capacity matrix
reallocv(RSTCKVEC(intordcm[ri][ci], gp));
reallocv(RSTCKVEC(intordcm[ri][ci], w));
gauss_points (gp.a, w.a, intordcm[ri][ci]);
for (i=0; i<intordcm[ri][ci]; i++)
{
if (ii == ipp)
{
xi = gp[i];
ncoord[0] = xi;
ncoord[1] = 0.0;
ncoord[2] = 0.0;
return 0;
}
ii++;
}
if (Tp->savemode==1) break;
}
if (Tp->savemode==1) break;
}
 
return 1;
}
 
 
 
/**
Function computes the total flux through a surface of an element
@param lcid - load case id
@param eid - element id (id in the list of all elements)
@param beid - id in the list of selected elements
@param flux - array of fluxes computed
TKr 06/08/2021 accroding to JK
*/
void linbart3d::surface_flux (long lcid,long eid,long beid,double *fluxes)
{
long i,ipp,fid;
double q,qn,area;
ivector nodes(nne);
bocontypet *bc;
vector flux(ASTCKVEC(ncomp)),a(ASTCKVEC(nne)),n(ASTCKVEC(ncomp));
// id of integration point for determination of the D matrix
ipp=Tt->elements[eid].ipp[0][0];
 
Tm->computenlfluxes (lcid,ipp);
Tm->givefluxes (lcid,ipp,flux);
 
// nodes on required elements
Tt->give_elemnodes (eid,nodes);
// array of cross section areas in nodes
Tc->give_area (eid,nodes,a);
 
 
// indicators of boundary conditions
bc = new bocontypet [nen];
Tb->bf[lcid].elemload[beid].give_bc (bc);
// loop over end nodes
for (i=0;i<nen;i++){
if (bc[i]==1){
// function constructs the outer normal vector
bar_normal_vectors (i,n);
// area of the i-th element surface
area = a[i];
// normal density flux
scprd (flux,n,qn);
// normal flux
q = qn*area;
// id of flux
fid = Tb->bf[lcid].elemload[beid].nvid[i];
// flux is stored to the appropriate position
fluxes[fid]+=q;
}
}
}
/trunk/SIFEL/TRFEL/SRC/linbart3d.h
0,0 → 1,114
#ifndef LINBART3D_H
#define LINBART3D_H
 
#include "genfile.h"
#include "aliast.h"
 
/**
class linbart3d defines onedimensional element with linear approximation functions
*/
class linbart3d
{
public:
linbart3d (void);
~linbart3d (void);
 
void dirvect (vector &s,vector &x,vector &y,vector &z);
void giveloccoord(vector &x,vector &y,vector &z,vector &lx);
 
void codnum (long *cn,long ri);
double approx (double xi,vector &nodval);
void intpointval (long eid);
void intpointval (long eid,vector &nodval,vector &ipval);
void initintpointval (long eid);
void intpointgrad (long eid);
void intpointother (long eid);
 
void btf_matrix (matrix &n,double xi,double v);
void bf_matrix (matrix &n,double xi);
void tf_matrix (matrix &w,double xi,double v);
void grad_matrix (matrix &gm,vector &x,vector &y,vector &z,double xi,double &jac);
void conductivity_matrix (long lcid,long eid,long ri,long ci,matrix &km);
void capacity_matrix (long eid,long ri,long ci,matrix &cm);
void quantity_source_vector (vector &sv,vector &nodval,long eid,long ri,long ci);
void transmission_matrix (long lcid,long eid,long ri,long ci,matrix &km);
void transmission_vector (vector &tmv,long lcid,long eid,long leid,long cid);
void convection_vector (vector &f,long lcid,long eid,long leid,long ri,long ci);
void internal_fluxes (long lcid,long eid,vector &ifl);
void accumul_energy (long lcid,long eid,vector &acum);
void advection_matrix (long lcid,long eid,long ri,long ci,matrix &hm);
 
void res_conductivity_matrix (long eid,long lcid,matrix &km);
void volume_rhs_vector (long lcid,long eid,long ri,long ci,vector &vrhs);
void volume_rhs_vector2 (long lcid,long eid,long ri,long ci,vector &vrhs);
void res_volume_rhs_vector (vector &f,long eid,long lcid);
void res_volume_rhs_vector2 (vector &f,long eid,long lcid);
void res_capacity_matrix (long eid,matrix &cm);
void res_convection_vector (vector &f,long lcid,long eid,long leid);
void res_transmission_vector (vector &f,long lcid,long eid,long leid);
void res_quantity_source_vector (vector &sv,vector &nodval,long lcid,long eid);
void res_internal_fluxes (long eid,vector &elemif);
void res_advection_matrix (long eid,long lcid,matrix &hm);
void nodal_energies (long lcid,long eid,vector &ifl);
void res_nodal_energy (long eid,vector &elemne,double dt);
/// function computes element quantity integral
/// the quantity is stored in nodes
double total_integral(long eid,vector &nodval);
/// function computes element quantity integral
/// the quantity is stored in integration points
double total_integral_ip (long eid,long varid);
 
 
void boundary_flux (vector &tmv,long lcid,long eid,long leid,long ri,long ci);
void res_boundary_flux (vector &f,long lcid,long eid,long leid);
void intpointflux (long eid);
void nod_grads_ip (long eid);
void nod_fluxes_ip (long eid);
void nod_others_comp (long lcid,long eid,long ri,long ci);
void transq_nodval (long eid,vector &nodval,nonmechquant nmq);
void transq_nodval_comp (long eid, vector &nodval, long ncne, long nq, nonmechquant *qt);
void transq_init_nodval (long eid,vector &nodval,nonmechquant nmq);
void transq_init_nodval_comp (long eid, vector &nodval, long ncne, long nq, nonmechquant *qt);
/// computes global coordinates of the given integration point of element
long ipcoord (long eid, long ipp, long ri, long ci, vector &coord);
/// function returns natural coordinates of the given integration point
long ipncoord (long eid, long ipp, vector &ncoord);
void surface_flux (long lcid,long eid,long beid,double *fluxes);
 
 
/// number of transported matters
long ntm;
/// total number of DOFs on the element
long ndofe;
/// numbers of DOFs for particular problems
long **dofe;
/// number of nodes on one element
long nne;
/// number of end nodes
long nen;
 
/// number of edges
long ned;
/// number of nodes on one edge
long nned;
 
/// number of approximated functions
long napfun;
/// problem dimension
long ncomp;
/// number of integration points
long **nip;
/// unknown ordering
long **ordering;
/// orders of integration of conductivity matrices
long **intordkm;
/// orders of integration of capacity matrices
long **intordcm;
};
 
#endif