BlueGreycalmElegant

Català-Valencià - Catalan中文 - Chinese (Simplified)中文 - Chinese (Traditional)Česky - CzechDansk - DanishNederlands - DutchEnglish - EnglishSuomi - FinnishFrançais - FrenchDeutsch - Germanעברית - Hebrewहिंदी - HindiMagyar - HungarianBahasa Indonesia - IndonesianItaliano - Italian日本語 - Japanese한국어 - KoreanМакедонски - Macedonianमराठी - MarathiNorsk - NorwegianPolski - PolishPortuguês - PortuguesePortuguês - Portuguese (Brazil)Русский - RussianSlovenčina - SlovakSlovenščina - SlovenianEspañol - SpanishSvenska - SwedishTürkçe - TurkishУкраїнська - UkrainianOëzbekcha - Uzbek

Compare Revisions

• This comparison shows the changes necessary to convert path

  → / @ 780

  → / @ 781

  ↔ Reverse comparison

• Compare Path:	Rev

  With Path:	Rev

Ignore whitespace Rev 780 → Rev 781

/trunk/SIFEL/GEFEL/gmatrix.cpp
1708,6 → 1708,10
scr -> mxv_scr(a,b);
break;
}
case symm_comp_columns:{
scc -> mxv_scc(a,b);
break;
}
case element_matrices:{
em -> mxv_em(a,b);
break;
1787,6 → 1791,10
scr -> addmat_scr (a,*gm.scr);
break;
}
case symm_comp_columns:{
scc -> addmat_scc (a,*gm.scc);
break;
}
case element_matrices:{
em -> addmat_em (a,*gm.em);
break;
1836,6 → 1844,10
scr -> scalmat_scr (a);
break;
}
case symm_comp_columns:{
scc -> scalmat_scc (a);
break;
}
case element_matrices:{
em -> scalmat_em (a);
break;

/trunk/SIFEL/MEFEL/PREP/hngen.cpp
0,0 → 1,13
#include <float.h>
#include "hngen.h"
#include "vector.h"
hngen::hngen()
{
x = y = z = 0.0;
xi = eta = zeta = 0.0;
cerr = DBL_MAX;
eid = -1;
et = gtypel(-1);
}

/trunk/SIFEL/MEFEL/PREP/hngen.h
0,0 → 1,20
#ifndef HNGEN_H
#define HNGEN_H
#include "galias.h"
#include "vector.h"
class hngen
{
public:
double x, y, z; // global coordinates of hanging node
double xi, eta, zeta; // natural coordinates of hanging node on element eid
double cerr; // error of computed natural coordinates if they are not in range [-1.0;1.0] exactly
long eid; // master element id on which the hanging node is located eid >=0,
gtypel et; // general element type used for the hanging node assignment
ivector mnodes; /// array of master nodes
hngen();
};
#endif

/trunk/SIFEL/MEFEL/SRC/dsm.cpp
0,0 → 1,1811
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "dsm.h"
#include "global.h"
#include "probdesc.h"
#include "mechmat.h"
#include "vecttens.h"
#include "intpoints.h"
#include "matrix.h"
#include "vector.h"
#include "tensor.h"
#include "elastisomat.h"
#include "mathem.h"
#include "nonlinman.h"
/**
This constructor inializes attributes to zero values.
*/
dsmmat::dsmmat (void)
{
m = 0.0;
lambda = 0.0;
kappa = 0.0;
c1_tilde = 0.0;
c2 = 0.0;
ks = 0.0;
//pr_suc_f = no;
//pr_tempr_f = no;
kappa_s0 = 0.0;
aks1 = 0.0;
aks2 = 0.0;
a0 = 0.0;
a2 = 0.0;
}
/**
This destructor is only for the formal purposes.
*/
dsmmat::~dsmmat (void)
{
}
/**
This function reads material parameters from the opened text file given
by the parameter in.
@param in - pointer to the opned text file
*/
void dsmmat::read (XFILE *in)
{
xfscanf (in, "%k%lf%k%lf%k%lf", "m", &m, "lambda", &lambda, "kappa", &kappa);
xfscanf (in, "%k%lf%k%lf%k%lf", "c1_tilde", &c1_tilde, "c2", &c2, "ks", &ks);
//preparation for the total stresses concept and the temperature effect:
//xfscanf (in, "%k%lf", "kappa_s0", &kappa_s0);
//xfscanf (in, "%k%lf%k%lf", "aks1", &aks1, "aks2", &aks2);
//xfscanf (in, "%k%lf%k%lf", "a0", &a0, "a2", &a2);
//xfscanf(in, "%k%m", "presc_suction", &answertype_kwdset, &pr_suc_f);
//if (pr_suc_fl == yes)
//suc_fn.read(in);
//xfscanf(in, "%k%m", "presc_tempr", &answertype_kwdset, &pr_tempr_f);
//if (pr_tempr_fl == yes)
//tempr_fn.read(in);
sra.read (in);
}
/**
This function initializes material model data
with respect of consistency parameter gamma.
@param ipp - integration point number
@param im - index of material type for given ip
@param ido - index of internal variables for given material in the ipp other array
12/06/2012 TKo
*/
void dsmmat::initval(long ipp, long im, long ido)
{
double v_ini, v_pc0, v_lambda1;
double i1s, j2s, s;
double v_kappa1, p1, bar_pc_0, pc_0;
long i, ncompstr = Mm->ip[ipp].ncompstr;
vector sig(ASTCKVEC(ncompstr)),sigt(ASTCKVEC(6)), q(ASTCKVEC(2));
double err=sra.give_err ();
// variables used for the calculation of initial value of preconsolidation pressure for the fully saturated state
double a, b, fs, ksi, sr, qq, p_atm = 101325.0;
if (Mm->ip[ipp].eqother[ido+ncompstr+1] == 0.0) // actual value of p_c was not set
{
// initial value of specific volume under small pressure p_1 on the normal consolidation line (v_lambda1 = 1+N)
v_lambda1 = Mm->ic[ipp][0];
// initial reference presure
p1 = Mm->ic[ipp][1];
// initial effective preconsolidation pressure \bar{p}_{c0}}
bar_pc_0 = Mm->ic[ipp][2];
// initial effective stresses
if ((Mp->eigstrains == 4) || (Mp->eigstrains == 5)){
for (i=0; i<ncompstr; i++)
sig(i) = Mm->ip[ipp].stress[i] = Mm->eigstresses[ipp][i];
}
else{
print_err("initial effective stresses (eigentsresses) are not defined on element %ld, ip=%ld,\n"
" cannot determine specific volume for initial stress state", __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1);
abort();
}
//vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
give_full_vector (sigt, sig, Mm->ip[ipp].ssst);
// initial value of mean stress
i1s = first_invar(sigt)/3.0;
q(0) = bar_pc_0;
s = Mm->givenonmechq(suction, ipp); // suction pressure s
sr = Mm->givenonmechq(saturation_degree, ipp); // degree of saturation S_r
//limitation of positive suction
if(s >=0.0)
s=0.0;
s = fabs(s);
q(1) = ks*s;
// test if the initial state is elastic, if not calculate new value of inital effective preconsolidation pressure pc_0 so the OCR=1
if (yieldfunction(sigt, q) > err){
j2s = j2_stress_invar (sigt);
bar_pc_0 = j2s/(m*m*(i1s-q(1))) + i1s;
}
//
// Calculation of inital value of preconsolidation pressure for the fully saturated state pc_0
//
// suction function:
//
// s / p_{atm}
// f(s) = 1 + -----------------------
// 10.7 + 2.4(s / p_{atm})
//
fs = 1.0 + (s/p_atm) / (10.7 + 2.4*(s/p_atm));
if (fs < 1.0)
fs = 1.0;
// bonding variable:
//
// \ksi = f(s) (1 - S_r)
//
ksi = (1-sr)*fs;
// according to Gallipoli
qq = 1.0 - c1_tilde*(1.0 - exp(c2*ksi)); // e/e_s
a = (qq - 1.0)/(lambda*qq - kappa)*(v_lambda1 - 1.0);
b = (lambda - kappa)/(lambda*qq - kappa);
pc_0 = -pow(-bar_pc_0 /exp(a), 1.0/b);
// compute initial value of specific volume
v_pc0 = v_lambda1 - lambda*log(pc_0/p1);
// actual value of p_c = initial effective preconsolidtaion pressure
Mm->ip[ipp].other[ido+ncompstr+1] = Mm->ip[ipp].eqother[ido+ncompstr+1] = bar_pc_0;
Mm->ip[ipp].other[ido+ncompstr+2] = Mm->ip[ipp].eqother[ido+ncompstr+2] = q(1);
// initial value of specific volume under small pressure p_1 on the swelling line (v_{\kappa l}) at the fully saturated state
Mm->ip[ipp].other[ido+ncompstr+3] = Mm->ip[ipp].eqother[ido+ncompstr+3] = v_kappa1 = v_pc0 + kappa*log(pc_0/p1);
// specific volume for intial stress state
Mm->ip[ipp].other[ido+ncompstr+4] = Mm->ip[ipp].eqother[ido+ncompstr+4] = v_ini = v_pc0 + kappa*log(pc_0/i1s);
// intial mean stress
Mm->ip[ipp].other[ido+ncompstr+5] = Mm->ip[ipp].eqother[ido+ncompstr+5] = i1s;
// saturation degree
Mm->ip[ipp].other[ido+ncompstr+9] = Mm->ip[ipp].eqother[ido+ncompstr+9] = Mm->givenonmechq(saturation_degree, ipp); // degree of saturation S_r
// initial depsv_dt
Mm->ip[ipp].other[ido+ncompstr+10] = Mm->ip[ipp].eqother[ido+ncompstr+10] = 0.0;
// porosity n for intial stress state
Mm->ip[ipp].other[ido+ncompstr+11] = Mm->ip[ipp].eqother[ido+ncompstr+11] = (v_ini-1.0)/v_ini; // v_ini = 1.0 + e_ini
// suction
Mm->ip[ipp].other[ido+ncompstr+12] = Mm->ip[ipp].eqother[ido+ncompstr+12] = Mm->givenonmechq(suction, ipp); // suction
}
// set inital value of Young's modulus
long idem = Mm->ip[ipp].gemid();
long ncompo = Mm->givencompeqother(ipp, im);
double e_ini = Mm->give_initial_ym(ipp, idem, ido+ncompo);
Mm->ip[ipp].other[ido+ncompstr+13] = Mm->ip[ipp].eqother[ido+ncompstr+13] = e_ini;
// time step scaling factor
Mm->ip[ipp].other[ido+ncompstr+14] = Mm->ip[ipp].eqother[ido+ncompstr+14] = 1.0;
// set inital stiffness matrix
matrix d;
makerefm(d, Mm->ip[ipp].other+ido+ncompstr+15, ncompstr, ncompstr);
Mm->elmatstiff(d, ipp, ido);
copym(d, Mm->ip[ipp].eqother+ido+ncompstr+15);
check_math_errel(Mm->elip[ipp]);
return;
}
/**
This function computes the value of yield functions.
@param sig - full stress tensor in Voigt notation (6 components)
@param q - %vector of hardening parameter
@retval The function returns value of yield function for the given stress tensor
25.3.2002
*/
double dsmmat::yieldfunction (vector &sig, vector &q)
{
double i1s,j2s,f;
i1s = first_invar (sig)/3.0;
j2s = j2_stress_invar (sig);
f = j2s/(m*m) + (i1s - q[1]) * (i1s - q[0]);
return f;
}
/**
This function computes derivatives of as-th yield function
with respect of vector sigma.
@param sig - stress tensor in Voigt notation (6 components)
@param q - %vector of the hardening parameter
@param dfds - %vector where the resulting derivatives are stored (6 components)
4.1.2002
*/
void dsmmat::deryieldfsigma (vector &sig, vector &q, vector &dfds)
{
double volume,i1s;
i1s = first_invar (sig)/3.0;
deviator (sig,dfds);
// q[0] = \bar{p}_c
// q[1] = p_s = k * s
cmulv(1.0/(m*m), dfds);
volume=1.0/3.0*(2.0*i1s - q[0] - q[1]);
dfds(0) += volume;
dfds(1) += volume;
dfds(2) += volume;
dfds(3) *= 2.0;
dfds(4) *= 2.0;
dfds(5) *= 2.0;
}
/**
The function computes the second derivatives of yield function
with respect of stress tensor sigma.
@param ddfds - tensor of the 4-th order where the are derivatives stored
19.12.2002
*/
void dsmmat::dderyieldfsigma (matrix &ddfds)
{
fillm(0.0, ddfds);
ddfds[0][0] = ddfds[1][1] = ddfds[2][2] = 2.0/(3.0*m*m) + 2.0/9.0;
ddfds[0][1] = ddfds[0][2] = ddfds[1][0] = ddfds[1][2] = -1.0/(3.0*m*m) + 2.0/9.0;
ddfds[2][0] = ddfds[2][1] = ddfds[0][1];
// doubled shear components due to tensor-> vector transformation
ddfds[3][3] = 2.0/(m*m);
ddfds[4][4] = 2.0/(m*m);
ddfds[5][5] = 2.0/(m*m);
}
/**
The function computes derivatives of plastic potential function
with respect of vector sigma.
@param sig - stress tensor
@param q - %vector of the hardening parameters
@param ddgds - %matrix where the resulting derivatives are stored
*/
void dsmmat::derpotsigma (vector &sig, vector &q, vector &ddgds)
{
deryieldfsigma (sig, q, ddgds);
}
/**
This function computes derivatives of as-th yield function
with respect of vector of hradening parameters.
@param sig - stress tensor
@param q - %vector of the hardening parameters
@param dfds - %vector where the resulting derivatives are stored
4.1.2002
*/
void dsmmat::deryieldfq(vector &sig, vector &q, vector &dfq)
{
// df
// ---------- = dfq[0]
// d\bar{p}_c
dfq[0] = -first_invar (sig)/3.0 + q[1];
// df
// --- = dfq[1]
// p_s
dfq[1] = -first_invar (sig)/3.0 + q[0];
return;
}
/**
The function computes the second derivatives of yield function
with respect to hradening parameters.
@param ddfdq - %matrix, where the resulting derivatives are stored
Created by Tomas Koudelka, 09.2012
*/
void dsmmat::deryieldfdqdq(matrix &ddfdq)
{
ddfdq[0][0] = 0.0;
ddfdq[0][1] = 1.0;
ddfdq[1][0] = 1.0;
ddfdq[1][1] = 0.0;
return;
}
/**
The function computes the second derivatives of yield function
with respect to stresses.
@param dfds - tensor, where the resulting derivatives are stored.
size of dfds = (6,number_of_hardening_param)
Created by Tomas Koudelka, 09.2012
*/
void dsmmat::deryieldfdsdq(matrix &dfdsdqt)
{
dfdsdqt[0][0] = dfdsdqt[1][0] = dfdsdqt[2][0] = -1.0/3.0;
dfdsdqt[0][1] = dfdsdqt[1][1] = dfdsdqt[2][1] = -1.0/3.0;
return;
}
/**
The function computes derivatives of hardening paramters
with respect of consistency parameter gamma.
@param ipp - integration point number
@param ido - index of internal variables for given material in the ipp other array
@param sig - full stress tensor
@param qtr - %vector of hardening variables
@param epsp - %vector of attained plastic strains
@param dqdg - %vector where the resulting derivatives are stored in Voigt notation
Created by Tomas Koudelka, 09.2012
*/
void dsmmat::der_q_gamma(long ipp, long ido, vector &sig, vector &qtr, vector &epsp, vector &dqdg)
{
double v_ini, v_lambda1, v_kappa1, p1, i1s, epsvp;
double p_atm = 101325.0;
double s, sr, fs, b, c, c1, ksi;
long ncompstr = Mm->ip[ipp].ncompstr;
strastrestate ssst = Mm->ip[ipp].ssst;
vector epst(ASTCKVEC(6));
// mean stress
i1s = first_invar (sig)/3.0;
// initial value of specific volume under small pressure p_1 for the normal consolidation line (v_{\lambda 1} = N)
v_lambda1 = Mm->ic[ipp][0];
// initial reference presure
p1 = Mm->ic[ipp][1];
// initial value of specific volume under small pressure p_1 for the swelling line (v_{\kappa 1})
v_kappa1 = Mm->ip[ipp].eqother[ido+ncompstr+3];
// specific volume for intial stress state
v_ini = Mm->ip[ipp].eqother[ido+ncompstr+4];
s = Mm->givenonmechq(suction, ipp); // suction pressure s
//limitation of positive suction
if(s >=0.0)
s=0.0;
s = fabs(s);
sr = Mm->givenonmechq(saturation_degree, ipp); // degree of saturation S_r
// actual value of preconsolidation pressure for the fully saturated state
give_full_vector (epst, epsp, ssst);
epsvp = first_invar(epst);
//double pc_new = p1 * exp((v_kappa1 - v_lambda1 + v_ini*epsvp)/(kappa-lambda));
// specific volume for the actual preconsolidation pressure:
//
// v_{pc} = v_{\lambda1} - \lambda \ln(p_c/p_1)
// double v_pc = v_lambda1 - lambda * log(pc_new/p1);
// suction function:
//
// s / p_{atm}
// f(s) = 1 + -----------------------
// 10.7 + 2.4(s / p_{atm})
//
fs = 1.0 + (s/p_atm) / (10.7 + 2.4*(s/p_atm));
if (fs < 1.0)
fs = 1.0;
// bonding variable:
//
// \ksi = f(s) (1 - S_r)
//
ksi = (1-sr)*fs;
// \tilde{c}_1
// c_1 = --------------------
// 1
// -------------- + 1
// v_{pc}-1
// c1 = c1_tilde/(1.0/(v_pc-1.0)+1.0); // according to Borja
c1 = c1_tilde; // according to Gallipoli
//
// c(\ksi) = 1 - c_1 [1-\exp(c_2 \ksi)]
//
c = 1.0 - c1*(1.0-exp(c2*ksi));
//
// \lambda - \kappa
// b(\ksi) = --------------------------
// \lambda c(\ksi) - \kappa
//
b = (lambda - kappa)/(lambda*c - kappa);
// q[0] = \bar{p}_c = -\exp[a(\ksi)] (-p_c)^{b(\ksi)}
//
// p_c = p_1 exp[(v_{\kappa_1} - v_{\lambda_1} + eps_{vp} v_{ini})/(\kappa - \lambda)]
//
// p_1 = const
// p_{c0} = const
// \sigma_{m0} = const
// v_{\kappa_1} = const
// v_{p_{c0}} = v_{\kappa_1} - \kappa \ln(p_{c0}/p_1) = const
// v_{ini} = v_{p_{c0}} + \kappa \ln(p_{c0}/\sigma_{m0}) = const
// v_{\lambda_1} = v_{p_{c0}} - \lambda \ln(p_{c0}/p_1) = const
//
// d \bar{p}_c d p_c (-p_c)^{b(\ksi)} d p_c
// ----------- = -\exp[a(\ksi)] b(\ksi)(-p_c)^{b(\ksi)-1} ---------- = -\exp[a(\ksi)] b(\ksi) ---------------- ---------- ,
// d \gamma d \gamma p_c d \gamma
//
// d p_c v_ini
// ---------- = p_c ----------------- (2p - \bar{p}_c - p_s),
// d \gamma \kappa - \lambda
//
// d \bar{p}_c b(\ksi) v_ini
// ----------- = -\bar{p}_c ---------------- (2p - \bar{p}_c - p_s)
// d \gamma \kappa - \lambda
//
//
dqdg[0] = -qtr[0] * (b*v_ini)/(kappa-lambda) * (2.0*i1s-qtr[0]-qtr[1]);
// q[1] = p_s = k*s
//
// d p_s
// ----------- = 0.0 because p_s is constant for the given strain level
// d \gamma
//
dqdg[1] = 0.0;
return;
}
/**
The function computes derivatives of hardening paramters
with respect of consistency parameter gamma.
@param ipp - integration point number
@param ido - index of internal variables for given material in the ipp other array
@param sig - full stress tensor
@param qtr - %vector of hardening variables
@param epsp - %vector of attained plastic strains
@param dqdgds - %matrix where the resulting derivatives are stored in Voigt notation,
its dimensions must be qtr.n x 6
Created by Tomas Koudelka, 05.2022
*/
void dsmmat::dqpardgammadsigma(long ipp, long ido, vector &sig, vector &qtr, vector &epsp, matrix &dqdgds)
{
double v_ini, v_lambda1, v_kappa1, p1, i1s, epsvp;
double p_atm = 101325.0;
double s, sr, fs, b, c, c1, ksi;
long ncompstr = Mm->ip[ipp].ncompstr;
strastrestate ssst = Mm->ip[ipp].ssst;
vector epst(ASTCKVEC(6));
// mean stress
i1s = first_invar (sig)/3.0;
// initial value of specific volume under small pressure p_1 for the normal consolidation line (v_{\lambda 1} = N)
v_lambda1 = Mm->ic[ipp][0];
// initial reference presure
p1 = Mm->ic[ipp][1];
// initial value of specific volume under small pressure p_1 for the swelling line (v_{\kappa 1})
v_kappa1 = Mm->ip[ipp].eqother[ido+ncompstr+3];
// specific volume for intial stress state
v_ini = Mm->ip[ipp].eqother[ido+ncompstr+4];
s = Mm->givenonmechq(suction, ipp); // suction pressure s
//limitation of positive suction
if(s >=0.0)
s=0.0;
s = fabs(s);
sr = Mm->givenonmechq(saturation_degree, ipp); // degree of saturation S_r
// actual value of preconsolidation pressure for the fully saturated state
give_full_vector (epst, epsp, ssst);
epsvp = first_invar(epst);
//double pc_new = p1 * exp((v_kappa1 - v_lambda1 + v_ini*epsvp)/(kappa-lambda));
// specific volume for the actual preconsolidation pressure:
//
// v_{pc} = v_{\lambda1} - \lambda \ln(p_c/p_1)
//double v_pc = v_lambda1 - lambda * log(pc_new/p1);
// suction function:
//
// s / p_{atm}
// f(s) = 1 + -----------------------
// 10.7 + 2.4(s / p_{atm})
//
fs = 1.0 + (s/p_atm) / (10.7 + 2.4*(s/p_atm));
if (fs < 1.0)
fs = 1.0;
// bonding variable:
//
// \ksi = f(s) (1 - S_r)
//
ksi = (1-sr)*fs;
// \tilde{c}_1
// c_1 = --------------------
// 1
// -------------- + 1
// v_{pc}-1
// c1 = c1_tilde/(1.0/(v_pc-1.0)+1.0); // according to Borja
c1 = c1_tilde; // according to Gallipoli
//
// c(\ksi) = 1 - c_1 [1-\exp(c_2 \ksi)]
//
c = 1.0 - c1*(1.0-exp(c2*ksi));
//
// \lambda - \kappa
// b(\ksi) = --------------------------
// \lambda c(\ksi) - \kappa
//
b = (lambda - kappa)/(lambda*c - kappa);
// q[0] = \bar{p}_c = -\exp[a(\ksi)] (-p_c)^{b(\ksi)}
//
// p_c = p_1 exp[(v_{\kappa_1} - v_{\lambda_1} + eps_{vp} v_{ini})/(\kappa - \lambda)]
//
// p_1 = const
// p_{c0} = const
// \sigma_{m0} = const
// v_{\kappa_1} = const
// v_{p_{c0}} = v_{\kappa_1} - \kappa \ln(p_{c0}/p_1) = const
// v_{ini} = v_{p_{c0}} + \kappa \ln(p_{c0}/\sigma_{m0}) = const
// v_{\lambda_1} = v_{p_{c0}} - \lambda \ln(p_{c0}/p_1) = const
//
// d \bar{p}_c d p_c (-p_c)^{b(\ksi)} d p_c
// ----------- = -\exp[a(\ksi)] b(\ksi)(-p_c)^{b(\ksi)-1} ---------- = -\exp[a(\ksi)] b(\ksi) ---------------- ---------- ,
// d \gamma d \gamma p_c d \gamma
//
// d p_c v_ini
// ---------- = p_c ----------------- (2p - \bar{p}_c - p_s),
// d \gamma \kappa - \lambda
//
// d \bar{p}_c b(\ksi) v_ini
// ----------- = -\bar{p}_c ---------------- (2p - \bar{p}_c - p_s)
// d \gamma \kappa - \lambda
//
//
// dqdg(0) = -qtr[0] * (b*v_ini)/(kappa-lambda) * (2.0*i1s-qtr[0]-qtr[1]);
dqdgds(0,0) = dqdgds(0,1) = dqdgds(0,2) = -qtr[0] * (b*v_ini)/(kappa-lambda) * 2.0/3.0;
// q[1] = p_s = k*s
//
// d p_s
// ----------- = 0.0 because p_s is constant for the given strain level
// d \gamma
//
dqdgds(1,0) = dqdgds(1,1) = dqdgds(1,2) = dqdgds(1,3) = dqdgds(1,4) = dqdgds(1,5) = 0.0;
return;
}
/**
The function computes derivatives of hardening paramters
with respect of consistency parameter gamma.
@param ipp - integration point number
@param ido - index of internal variables for given material in the ipp other array
@param sig - full stress tensor
@param qtr - %vector of hardening variables
@param epsp - %vector of attained plastic strains
@param dqdgds - %matrix where the resulting derivatives are stored, its dimensions must be qtr.n x qtr.n
Created by Tomas Koudelka, 05.2022
*/
void dsmmat::dqpardgammadqpar(long ipp, long ido, vector &sig, vector &qtr, vector &epsp, matrix &dqdgdq)
{
double v_ini, v_lambda1, v_kappa1, p1, i1s, epsvp;
double p_atm = 101325.0;
double s, sr, fs, b, c, c1, ksi;
long ncompstr = Mm->ip[ipp].ncompstr;
strastrestate ssst = Mm->ip[ipp].ssst;
vector epst(ASTCKVEC(6));
// mean stress
i1s = first_invar (sig)/3.0;
// initial value of specific volume under small pressure p_1 for the normal consolidation line (v_{\lambda 1} = N)
v_lambda1 = Mm->ic[ipp][0];
// initial reference presure
p1 = Mm->ic[ipp][1];
// initial value of specific volume under small pressure p_1 for the swelling line (v_{\kappa 1})
v_kappa1 = Mm->ip[ipp].eqother[ido+ncompstr+3];
// specific volume for intial stress state
v_ini = Mm->ip[ipp].eqother[ido+ncompstr+4];
s = Mm->givenonmechq(suction, ipp); // suction pressure s
//limitation of positive suction
if(s >=0.0)
s=0.0;
s = fabs(s);
sr = Mm->givenonmechq(saturation_degree, ipp); // degree of saturation S_r
// actual value of preconsolidation pressure for the fully saturated state
give_full_vector (epst, epsp, ssst);
epsvp = first_invar(epst);
//double pc_new = p1 * exp((v_kappa1 - v_lambda1 + v_ini*epsvp)/(kappa-lambda));
// specific volume for the actual preconsolidation pressure:
//
// v_{pc} = v_{\lambda1} - \lambda \ln(p_c/p_1)
//double v_pc = v_lambda1 - lambda * log(pc_new/p1);
// suction function:
//
// s / p_{atm}
// f(s) = 1 + -----------------------
// 10.7 + 2.4(s / p_{atm})
//
fs = 1.0 + (s/p_atm) / (10.7 + 2.4*(s/p_atm));
if (fs < 1.0)
fs = 1.0;
// bonding variable:
//
// \ksi = f(s) (1 - S_r)
//
ksi = (1-sr)*fs;
// \tilde{c}_1
// c_1 = --------------------
// 1
// -------------- + 1
// v_{pc}-1
// c1 = c1_tilde/(1.0/(v_pc-1.0)+1.0); // according to Borja
c1 = c1_tilde; // according to Gallipoli
//
// c(\ksi) = 1 - c_1 [1-\exp(c_2 \ksi)]
//
c = 1.0 - c1*(1.0-exp(c2*ksi));
//
// \lambda - \kappa
// b(\ksi) = --------------------------
// \lambda c(\ksi) - \kappa
//
b = (lambda - kappa)/(lambda*c - kappa);
// q[0] = \bar{p}_c = -\exp[a(\ksi)] (-p_c)^{b(\ksi)}
//
// p_c = p_1 exp[(v_{\kappa_1} - v_{\lambda_1} + eps_{vp} v_{ini})/(\kappa - \lambda)]
//
// p_1 = const
// p_{c0} = const
// \sigma_{m0} = const
// v_{\kappa_1} = const
// v_{p_{c0}} = v_{\kappa_1} - \kappa \ln(p_{c0}/p_1) = const
// v_{ini} = v_{p_{c0}} + \kappa \ln(p_{c0}/\sigma_{m0}) = const
// v_{\lambda_1} = v_{p_{c0}} - \lambda \ln(p_{c0}/p_1) = const
//
// d \bar{p}_c d p_c (-p_c)^{b(\ksi)} d p_c
// ----------- = -\exp[a(\ksi)] b(\ksi)(-p_c)^{b(\ksi)-1} ---------- = -\exp[a(\ksi)] b(\ksi) ---------------- ---------- ,
// d \gamma d \gamma p_c d \gamma
//
// d p_c v_ini
// ---------- = p_c ----------------- (2p - \bar{p}_c - p_s),
// d \gamma \kappa - \lambda
//
// d \bar{p}_c b(\ksi) v_ini
// ----------- = -\bar{p}_c ---------------- (2p - \bar{p}_c - p_s)
// d \gamma \kappa - \lambda
//
//
// dqdg(0) = -qtr[0] * (b*v_ini)/(kappa-lambda) * (2.0*i1s-qtr[0]-qtr[1]);
// wrong //dqdgdq(0,0) = 2.0*qtr[0] * (b*v_ini)/(kappa-lambda);
dqdgdq(0,0) = -(b*v_ini)/(kappa-lambda) * (2.0*i1s-2.0*qtr[0]-qtr[1]);
dqdgdq(0,1) = qtr[0] * (b*v_ini)/(kappa-lambda);
// q[1] = p_s = k*s
//
// d p_s
// ----------- = 0.0 because p_s is constant for the given strain level
// d \gamma
//
dqdgdq(1,0) = dqdgdq(1,1) = 0.0;
return;
}
/**
This function computes plastic modulus.
@param ipp - integration point number
@param ido - index of internal variables for given material in the ipp other array
@param sig - stress tensor in Voigt notation
@param qtr - %vector of hardening parameters
@param epsp - %vector of attained plastic strains
@retval The function returns value of the plastic modulus.
Created by Tomas Koudelka, 09.2012
*/
double dsmmat::plasmodscalar(long ipp, long ido, vector &sig, vector &qtr, vector &epsp)
{
double ret;
vector dfq(ASTCKVEC(2));
vector dqg(ASTCKVEC(2));
deryieldfq(sig, qtr, dfq);
der_q_gamma(ipp, ido, sig, qtr, epsp, dqg);
scprd(dfq, dqg, ret);
return -ret;
}
/**
This function computes new value of the hardening parameter q.
@param ipp - integration point pointer
@param ido - index of internal variables for given material in the ipp other array
@param eps - %vector of the attained strains in Voigt notation
@param epsp - %vector of the attained plastic strains in Voigt notation
@param sig - attained stress tensor in Voigt notation
@param ssst - stress/strain state parameter (for the used vector_tensor function)
@param q - %vector of the hardening parameters
Created by Tomas Koudelka, 09.2012
*/
void dsmmat::updateq(long ipp, long ido, vector &eps, vector &epsp, vector &sig, vector &q)
{
vector epst(ASTCKVEC(6));
vector depsp(epsp.n);
double v_kappa1, v_lambda1, p1, depsvp, v_ini, pc_new;
double v_pc, p_atm = 101325.0;
double s, sr, fs, a, b, c, c1, ksi;
strastrestate ssst = Mm->ip[ipp].ssst;
long ncompstr = Mm->ip[ipp].ncompstr;
long i;
// initial value of specific volume under small pressure p_1 for the normal consolidation line (v_lambda1 = 1+N)
v_lambda1 = Mm->ic[ipp][0];
// initial reference presure
p1 = Mm->ic[ipp][1];
// original specific volume before any loading
v_kappa1 = Mm->ip[ipp].eqother[ido+ncompstr+3];
// specific volume for intial stress state
v_ini = Mm->ip[ipp].eqother[ido+ncompstr+4];
// Original version of pc calculation - problems with tension
// vk = v_ini*(1+epsv)+kappa*log(i1s/p1);
// pc_new = p1 * exp((vk-v_lambda1)/(kappa-lambda));
//
// Alternative version of pc calculation - 'less' accurate in benchmarks
// pc_new = p1 * exp((v_kappa1 - v_lambda1 + v_ini*epsv)/(kappa-lambda));
// vk = v_ini*(1+epsv)+kappa*log(i1s/p1);
for (i=0; i<ncompstr; i++)
depsp[i] = epsp[i];// + Mm->ic[ipp][3+i];
give_full_vector(epst, depsp, ssst);
depsvp = first_invar(epst);
pc_new = p1 * exp((v_kappa1 - v_lambda1 + v_ini*depsvp)/(kappa-lambda));
// For the partially saturated soil,
// the effective preconsolidation pressure \bar{p}_c is related to the saturated
// preconsolidation pressure p_c according to bonding variable \ksi
//
// \bar{p}_c = -exp[a(\ksi)] (-p_c)^{b(\ksi)}
//
//
// specific volume for the actual preconsolidation pressure:
//
// v_{pc} = v_{\lambda1} - \lambda \ln(p_c/p_1)
v_pc = v_lambda1 - lambda * log(pc_new/p1);
s = Mm->givenonmechq(suction, ipp); // suction pressure s
//limitation of positive suction
if(s >=0.0)
s=0.0;
s = fabs(s);
sr = Mm->givenonmechq(saturation_degree, ipp); // degree of saturation S_r
// suction function:
//
// s / p_{atm}
// f(s) = 1 + -----------------------
// 10.7 + 2.4(s / p_{atm})
//
fs = 1.0 + (s/p_atm) / (10.7 + 2.4*(s/p_atm));
if (fs < 1.0)
fs = 1.0;
// bonding variable:
//
// \ksi = f(s) (1 - S_r)
//
ksi = (1-sr)*fs;
// \tilde{c}_1
// c_1 = --------------------
// 1
// -------------- + 1
// v_{pc}-1
// c1 = c1_tilde/(1.0/(v_pc-1.0)+1.0); // according to Borja
c1 = c1_tilde; // according to Gallipoli
//
// c(\ksi) = 1 - c_1 [1-\exp(c_2 \ksi)]
//
c = 1.0 - c1*(1.0-exp(c2*ksi));
//
// a(\ksi) = v_{lambda_1} [c(\ksi)-1]/[\lambda c(\ksi) - \kappa)
//
a = (v_lambda1 - 1.0)*(c - 1.0)/(lambda*c - kappa);
//
// \lambda - \kappa
// b(\ksi) = --------------------------
// \lambda c(\ksi) - \kappa
//
b = (lambda - kappa)/(lambda*c - kappa);
//
// effective preconsolidation pressure:
//
// q[0] = \bar{p}_c = -exp[a(\ksi)] (-p_c)^{b(\ksi)}
//
q[0] = -exp(a)*pow(-pc_new, b);
// apparent adhesion due to suction
// (right intersection of the yield surface with p axis):
//
// q[1] = p_s = k s
//
q[1] = ks * s;
return;
}
/**
This function computes material stiffnes matrix.
@param d - allocated matrix structure for material stiffness %matrix in the reduced format
@param ipp - integration point number
@param ido - index of internal variables for given material in the ipp other array
Created by Tomas Koudelka, 09.2012
*/
void dsmmat::matstiff (matrix &d, long ipp, long ido)
{
double zero=1.0e-20;
switch (Mp->nlman->stmat){
case initial_stiff:
// initial elastic matrix
Mm->elmatstiff(d, ipp, ido);
break;
case tangent_stiff:
case secant_stiff:
case incr_tangent_stiff:
case ijth_tangent_stiff:{
// consistent tangent stiffness matrix
long n = Mm->ip[ipp].ncompstr;
if (sra.give_tsra() == cp){
// algorithmic consistent stiffness matrix for the cutting plane algorithm is computed in the course of stress return
if (d.m != n){ // material stiffness matrix used for element in the case of plane stress/strain states
for(long i=0; i<d.m; i++) // copy just the first d.n(=3) components at all matrix rows
memcpy(&d(i,0), Mm->ip[ipp].other+ido+n+15+n*i, d.n*sizeof(*Mm->ip[ipp].other));
}
else{ // dimensions of material stiffness matrix matches the ones for element material stiffness matrix
matrix auxm;
makerefm(auxm, Mm->ip[ipp].other+ido+n+15, n, n);
copym(auxm, d);
}
}
else{
// compute algorithmic stiffness matrix for closest point projection stress return algorithm
matrix de(ASTCKMAT(6, 6)), ddfdst(ASTCKMAT(6,6));
matrix e(ASTCKMAT(6, 6)), auxd(ASTCKMAT(6,6)), tmp(ASTCKMAT(6, 6));
double gamma = Mm->ip[ipp].other[ido+n];
identm(e);
Mm->elmatstiff(de, ipp, spacestress);
dderyieldfsigma(ddfdst);
cmulm(gamma, de, auxd);
mxm(auxd, ddfdst, tmp);
addm(e, tmp, auxd);
invm(auxd, tmp, zero);
mxm(tmp, de, auxd);
// matrix representation of the fourth order tensor
tensor4_ematrix(d, auxd, Mm->ip[ipp].ssst);
}
break;
}
}
}
/**
This function computes elasto-plastic material stiffnes matrix.
@param d - allocated matrix structure for material stiffness %matrix in the reduced format
@param ipp - integration point number
@param ido - index of internal variables for given material in the ipp other array
Created by Tomas Koudelka, 09.2012
*/
void dsmmat::depmatstiff (matrix &d, long ipp, long ido)
{
switch (Mp->nlman->stmat){
case initial_stiff:
// initial elastic matrix
Mm->elmatstiff(d, ipp, ido);
break;
case tangent_stiff:
case secant_stiff:
case incr_tangent_stiff:
case ijth_tangent_stiff:{
// consitent tangent stiffness matrix
long n = Mm->ip[ipp].ncompstr;
long ncompq = 2;
matrix de(ASTCKMAT(n, n)), auxm;
vector sig(ASTCKVEC(n)), q(ASTCKVEC(ncompq)), epsp(ASTCKVEC(n));
vector dfds(ASTCKVEC(n)), dgds(ASTCKVEC(n));
vector sigt(ASTCKVEC(6)), dfdst(ASTCKVEC(6)), dgdst(ASTCKVEC(6));
vector dedgds(ASTCKVEC(n)), dedfds(ASTCKVEC(n));
strastrestate ssst = Mm->ip[ipp].ssst;
double denomh, denomd, denom, gamma, f, err=sra.give_err();
// prepare actual values of needed quantities
copyv(Mm->ip[ipp].stress, sig); // stresses
copyv(Mm->ip[ipp].other+ido, epsp); // plastic strains
gamma = Mm->ip[ipp].other[ido+n+0]; // cumulated plastic multiplier
q(0) = Mm->ip[ipp].other[ido+n+1]; // hardening parameter \bar{p}_c
q(1) = Mm->ip[ipp].other[ido+n+2]; // hardening parameter k_s
give_full_vector(sigt, sig, ssst); // 6 component vector of stresses
// check actual stress state: f <= 0.0
f = yieldfunction(sigt, q);
if ((f < -err) || (gamma == 0.0 && (f < err))){
// stress state is in the elastic domain,
// i.e. material has been subjected by elastic loading/unloading ->
// -> return actual elastic stiffness matrix
Mm->elmatstiff(d, ipp, ssst);
return;
}
// stress state is on/close to the yield surface
// material has been subjected by plastic loading ->
// compute elasto-plastic matrix
// D_e \frac{dg}{d\sigma} (\frac{df}{d\sigma})^T D_e
// D_{ep} = D_e - ---------------------------------------------------
// (\frac{df}{d\sigma})^T D_e \frac{dg}{d\sigma} + H
//
// H = -\frac{df}{dq} \frac{dq}{d\lambda}
Mm->elmatstiff(de, ipp, ssst);
deryieldfsigma(sigt, q, dfdst);
derpotsigma(sigt, q, dgdst);
// conversion of full tensor notation to vector notation
give_red_vector(dfdst, dfds, ssst);
// conversion of full tensor notation to vector notation
give_red_vector(dgdst, dgds, ssst);
mxv(de,dgds,dedgds);
mxv(de,dfds,dedfds);
denomh = plasmodscalar(ipp, ido, sig, q, epsp);
scprd(dfds, dedgds, denomd);
// preserve minimum stiffness for perfect plastic material, i.e. if the hardening modulus is 'close' to zero (H->0)
// if (fabs(denomh/denomd) < 1.0e-3)
// denomh = sgn(denomh)*1.0e-3*denomd;
denom = denomd+denomh;
if (n == d.m){
// material stiffness matrix has identical dimensions with the ones of stiffness matrix used in
// assembling of element stiffness matrix on finite element -> auxm will be reference on matrix d
makerefm(auxm, d.a, d.m, d.n);
nullm(auxm);
}
else{
// material stiffness matrix has different dimensions than the stiffness matrix used in
// assembling of element stiffness matrix on finite element -> auxm will be reference on matrix d
reallocm(RSTCKMAT(n, n, auxm));
}
// calculate matrix in the nominator
vxv(dedgds, dedfds, auxm);
// scaling by the denominator
cmulm(1.0/denom, auxm);
// resulting material stiffness matrix
subm(de, auxm, auxm);
// resulting material stiffness matrix with dimensions adjusted for finite element computation
if (n != d.m){
// resulting material stiffness matrix dimensions must be adjusted for finite element computation
extractblock(auxm, d, 0, 0);
}
else{
// nothing is needed, matrix auxm is direct reference to the matrix in argument d
// hence d contains resulting material stiffness matrix with appropriate dimensions for
// the stiffness matrix on the finite element
}
/*
double zero=1.0e-20;
matrix de(ASTCKMAT(6, 6)), ddfdst(ASTCKMAT(6,6));
matrix e(ASTCKMAT(6, 6)), auxd(ASTCKMAT(6,6)), tmp(ASTCKMAT(6, 6));
long n = Mm->ip[ipp].ncompstr;
double gamma = Mm->ip[ipp].other[ido+n];
identm(e);
Mm->elmatstiff(de, ipp, spacestress);
dderyieldfsigma(ddfdst);
cmulm(gamma, de, auxd);
mxm(auxd, ddfdst, tmp);
addm(e, tmp, auxd);
invm(auxd, tmp, zero);
mxm(tmp, de, auxd);
// matrix representation of the fourth order tensor
tensor4_ematrix(d, auxd, Mm->ip[ipp].ssst);*/
break;
}
}
}
/**
The function computes stress increments due to suction and temperature changes in the integration point and stores
them into eqother array.
@param ipp - integration point pointer
@param im - index of material type for given ip
@param ido - index of internal variables for given material in the ipp other array
@return The function does not return anything but stores the resulting values in the eqother array.
Created by Tomas krejci, 03/2023
*/
void dsmmat::nlstressesincr (long ipp, long im, long ido)
//
{
double n,v,depss=0.0,depst=0.0,ds,s,kappa_s,lambda_s,p_atm = 101325.0,p,pc_ref;
double temp=0.0,temp0=0.0,dtemp=0.0;
long i,ncomp;
// stress/strain state
strastrestate ss = Mm->ip[ipp].ssst;
ncomp=Mm->ip[ipp].ncompstr;
matrix d(ASTCKMAT(ncomp,ncomp));
vector depsv(ASTCKVEC(ncomp)),dsigma(ASTCKVEC(ncomp));
//if (Mm->ip[ipp].tm[im-1] == effective_stress){
// In case of effective stresses concept, stress increments due to suction changes are included in eff_stress material
//for (i=0;i<ncomp;i++)
Mm->ip[ipp].eqother[ido+ncomp+15+ncomp*ncomp+i] = 0.0;
//}
//else{
//}
}
/**
The function computes stress increment due to pore pressure change in the integration point and stores
them into ip stress array.
@param lcid - load case id
@param ipp - integration point number
@param im - index of the material in the tm and idm arrays on integration point
@param ido - index of internal variables in the ip's ipp eqother array
@param fi - first index of the required stress increment component
@param sig - %vector containing stress increment components (output)
@return The function returns %vector of stress increments in the parameter sig.
Created by Tomas Krejci, 03/2023
*/
void dsmmat::givestressincr (long lcid, long ipp, long im, long ido, long fi, vector &sig)
{
nullv(sig);
long i, j, ncomp = Mm->ip[ipp].ncompstr;
for (i=fi,j=0; j<sig.n; i++, j++)
sig(j) = Mm->ip[ipp].eqother[ido+ncomp+15+ncomp*ncomp+i];
}
/**
This function computes stresses at given integration point ipp,
depending on the reached strains.
The cutting plane algorithm is used. The stress and the other attribute of
given integration point is actualized.
@param ipp - integration point number in the mechmat ip array.
@param im - index of material type for given ip
@param ido - index of internal variables for given material in the ipp other array
Created by Tomas Koudelka, 09.2012
*/
void dsmmat::nlstresses (long ipp, long im, long ido)
//
{
int errr;
long i,ni,ncomp=Mm->ip[ipp].ncompstr, matcomp=0, ret;
double gamma, err, pnewdt = 1.0;
vector epsn(ASTCKVEC(ncomp)), epsp(ASTCKVEC(ncomp)), q(ASTCKVEC(2));
vector epsa(ASTCKVEC(ncomp)), sig(ASTCKVEC(ncomp)), dfds(ASTCKVEC(ncomp)), dgds(ASTCKVEC(ncomp));
matrix d;
vector dfdst(ASTCKVEC(6)), dgdst(ASTCKVEC(6));
vector epst(ASTCKVEC(6)), epsnt(ASTCKVEC(6)), epspt(ASTCKVEC(6)), sigt(ASTCKVEC(6));
//vector dev(ASTCKVEC(6));
double i1s, j2s, epsv, epsvp;
double dtime = Mp->timecon.actualforwtimeincr();
// stress/strain state
strastrestate ssst = Mm->ip[ipp].ssst;
double nu;
// initial values
for (i=0; i<ncomp; i++)
{
epsn[i] = Mm->ip[ipp].strain[i];
epsp[i] = Mm->ip[ipp].eqother[ido+i];
// It is necessary to substract initial strains in order to invoke initial nonzero stress
// state for zero displacements (i.e. zero strains)
// sig = E*(eps - eps_p), sig_0 = E*eps_0. For eps = 0: sig = sig_0 => eps_p = -eps_0
// epsp[i] -= Mm->ic[ipp][3+i];
}
gamma=Mm->ip[ipp].eqother[ido+ncomp];
q[0] = Mm->ip[ipp].eqother[ido+ncomp+1]; // effective preconsolidation pressure \bar{p}_c
q[1] = Mm->ip[ipp].eqother[ido+ncomp+2]; // apparent adhesion p_s = k * s
// stress return algorithm
switch(sra.give_tsra ())
{
case cp:
ni=sra.give_ni ();
err=sra.give_err ();
matcomp = actualize_stiffmat(*Mp->nlman, Mp->istep, Mp->jstep);
if (Mp->istep >= 265){
matcomp = 0;
}
if (matcomp) reallocm(RSTCKMAT(ncomp, ncomp, d));
ret = Mm->cutting_plane (ipp, im, ido, gamma, epsn, epsp, q, ni, err, matcomp, d);
if (ret)
pnewdt = 0.5;
break;
case gsra:
ni=sra.give_ni ();
err=sra.give_err ();
Mm->newton_stress_return (ipp, im, ido, gamma, epsn, epsp, q, ni, err);
break;
default:
print_err("wrong type of stress return algorithm is required", __FILE__, __LINE__, __func__);
abort ();
}
/*
////////////////////////////////
// elastic material for debug??
// initial values
for (i=0;i<ncomp;i++){
epsn[i]=epsp[i]=Mm->ip[ipp].strain[i];
}
reallocm(RSTCKMAT(ncomp, ncomp, d));
Mm->elmatstiff (d,ipp,ssst);
nu = 0.3;
mxv (d,epsn,sig);
if (Mm->ip[ipp].ssst == planestress)
Mm->ip[ipp].strain[3] = -nu / (1.0 - nu) * (epsn[0]+epsn[1]);
if ((Mp->eigstrains == 4) || (Mp->eigstrains == 5))
{
for (i=0;i<ncomp;i++)
sig(i) += Mm->eigstresses[ipp][i];
}
for (i=0;i<ncomp;i++){
Mm->ip[ipp].stress[i]=sig(i);
}
////////////////////////////////
*/
// new data storage
for (i=0; i<ncomp; i++){
//epsp[i] += Mm->ic[ipp][3+i];
Mm->ip[ipp].other[ido+i]=epsp[i];
}
Mm->ip[ipp].other[ido+ncomp]=gamma;
Mm->ip[ipp].other[ido+ncomp+1]=q[0]; // the first hardening parameter \bar{p}_c
Mm->ip[ipp].other[ido+ncomp+2]=q[1]; // the second hardening parameter p_s
// other[ido+ncomp+3] = v_lambda1; it is the initial value and therefore it does not require actualization
// other[ido+ncomp+4] = v_ini; it is the initial value and therefore it does not require actualization
Mm->givestress(0, ipp, sig);
give_full_vector(sigt,sig,Mm->ip[ipp].ssst);
give_full_vector(epsnt,epsn,Mm->ip[ipp].ssst);
give_full_vector(epspt,epsp,Mm->ip[ipp].ssst);
i1s = first_invar (sigt)/3.0;
// the second invariant of deviator
// it is expressed with the help of the stress components
// components of the deviator are not needed
j2s = sqrt(j2_stress_invar (sigt));
epsv = first_invar (epsnt);
epsvp = first_invar (epspt);
Mm->ip[ipp].other[ido+ncomp+5] = i1s;
Mm->ip[ipp].other[ido+ncomp+6] = j2s;
Mm->ip[ipp].other[ido+ncomp+7] = epsv;
Mm->ip[ipp].other[ido+ncomp+8] = epsvp;
Mm->ip[ipp].other[ido+ncomp+9] = Mm->givenonmechq(saturation_degree, ipp);
// volumetric strain rate
//Mm->ip[ipp].other[ido+ncomp+10] = (epsv - Mm->ip[ipp].eqother[ido+ncomp+7])/dtime;
// the volumetric strain rate is computed via generalized trapesoidal rule
Mm->ip[ipp].other[ido+ncomp+10] = 0.5*((epsv - Mm->ip[ipp].eqother[ido+ncomp+7])/dtime + Mm->ip[ipp].eqother[ido+ncomp+10]);
// reference presure = p1
double p1 = Mm->ic[ipp][1];
// specific volume at the reference pressure p_1 = v_lambda1
double v_lambda1 = Mm->ic[ipp][0];
// specific volume for the actual preconsolidation pressure = v_pc
double v_pc = v_lambda1 - lambda * log(q[0]/p1);
// specific volume for the actual pressure p
double v;
if (i1s < 0.0)
v = v_pc + kappa*log(q[0]/i1s);
else
v = v_pc + kappa*log(-q[0]);
// actual porosity n
Mm->ip[ipp].other[ido+ncomp+11] = (v-1.0)/v;
// suction
Mm->ip[ipp].other[ido+ncomp+12] = Mm->givenonmechq(suction, ipp);
// Young modulus
Mm->ip[ipp].other[ido+ncomp+13] = give_actual_ym(ipp, im, ido);
// time step scaling factor
Mm->ip[ipp].other[ido+ncomp+14] = pnewdt;
if (matcomp)
copym(d, Mm->ip[ipp].other+ido+ncomp+15);
errr = check_math_errel(Mm->elip[ipp]);
}
/**
Function returns flag whether actualize the stiffness matrix according to the setup of the nonlinear solver and
attained number of steps.
@param[in] smt - type of the stiffness matrix modification in NR method
@param[in] istep - number of load increment step in NR method
@param[in] jstep - number of inner iteration step in NR method
@retval 0 - in the case that actualization of the stiffness matrix is NOT required
@retval 1 - in the case that actualization of the stiffness matrix is required
*/
long dsmmat::actualize_stiffmat(const nonlinman &nlman, long istep, long jstep)
{
long ret = 0;
switch (nlman.stmat){
case initial_stiff:
break;
case tangent_stiff:
ret = 1;
break;
case secant_stiff:
case incr_tangent_stiff:{
if (jstep <= 0)
ret = 1;
break;
}
case ijth_tangent_stiff:{
if (((jstep <= 0) && (istep % nlman.ithstiffchange == 0)) ||
((jstep+1) % nlman.jthstiffchange == 0))
ret = 1;
break;
}
default:{
print_err("unknown type of stiffness matrix is required",__FILE__,__LINE__,__func__);
}
}
return ret;
}
/**
This function updates values in the other array reached in the previous equlibrium state to
values reached in the new actual equilibrium state.
@param ipp - integration point number in the mechmat ip array.
@param im - index of material type for given ip
@param ido - index of internal variables for given material in the ipp other array
Created by Tomas Koudelka, 09.2012
*/
void dsmmat::updateval (long ipp, long im, long ido)
{
long i,n = Mm->ip[ipp].ncompstr;
double e = comp_actual_ym(ipp, im, ido);
for (i=0;i<n;i++){
Mm->ip[ipp].eqother[ido+i]=Mm->ip[ipp].other[ido+i];
}
Mm->ip[ipp].eqother[ido+n] = Mm->ip[ipp].other[ido+n]; // gamma
Mm->ip[ipp].eqother[ido+n+1] = Mm->ip[ipp].other[ido+n+1]; // p_c (hardening parameter - effective preconsolidation pressure)
Mm->ip[ipp].eqother[ido+n+2] = Mm->ip[ipp].other[ido+n+2]; // p_s (hardening parameter - apparent adhesion)
//
// Following three values of the specific volume are constant and they
// are stored directly in the eqother array during the intial step
//
// Mm->ip[ipp].eqother[ido+n+3] = Mm->ip[ipp].other[ido+n+3]; // v_lambda1
// Mm->ip[ipp].eqother[ido+n+4] = Mm->ip[ipp].other[ido+n+4]; // v_ini
//
Mm->ip[ipp].eqother[ido+n+5] = Mm->ip[ipp].other[ido+n+5]; // i1s;
Mm->ip[ipp].eqother[ido+n+6] = Mm->ip[ipp].other[ido+n+6]; // j2s;
Mm->ip[ipp].eqother[ido+n+7] = Mm->ip[ipp].other[ido+n+7]; // epsv;
Mm->ip[ipp].eqother[ido+n+8] = Mm->ip[ipp].other[ido+n+8]; // epsvp;
// saturation degree
Mm->ip[ipp].eqother[ido+n+9] = Mm->ip[ipp].other[ido+n+9]; // Sr
// volumteric strain rate
Mm->ip[ipp].eqother[ido+n+10] = Mm->ip[ipp].other[ido+n+10]; // depsv/dt
// actual porosity n
Mm->ip[ipp].eqother[ido+n+11] = Mm->ip[ipp].other[ido+n+11];
// suction
Mm->ip[ipp].eqother[ido+n+12] = Mm->ip[ipp].other[ido+n+12]; //s
// Young modulus
Mm->ip[ipp].eqother[ido+n+13] = Mm->ip[ipp].other[ido+n+13] = e;
// time step scaling factor
Mm->ip[ipp].eqother[ido+n+14] = Mm->ip[ipp].other[ido+n+14];
// stiffness matrix
for (i=0;i<n*n;i++){
Mm->ip[ipp].eqother[ido+n+15+i] = Mm->ip[ipp].other[ido+n+15+i];
}
// stress increments - it is not needed
//for (i=0;i<n;i++){
//Mm->ip[ipp].eqother[ido+n+15+n*n+i]=Mm->ip[ipp].other[ido+n+15+n*n+i];
//}
}
/**
The function returns actual Young modulus value.
@param ipp - integration point number
@param im - index of material type for given ip
@param ido - index of internal variables for given material in the ipp other array
*/
double dsmmat::give_actual_ym(long ipp, long im, long ido)
{
long ncomp=Mm->ip[ipp].ncompstr;
double e = Mm->ip[ipp].other[ido+ncomp+13];
return e;
}
/**
The function computes actual Young modulus value.
@param ipp - integration point number
@param im - index of material type for given ip
@param ido - index of internal variables for given material in the ipp other array
*/
double dsmmat::comp_actual_ym(long ipp, long im, long ido)
{
long idem = Mm->ip[ipp].gemid();
long ncompo = Mm->givencompeqother(ipp, im);
double e;
// initial Young modulus
double e_ini = Mm->give_initial_ym(ipp, idem, ido+ncompo);
// actual Poisson's ratio
double nu = give_actual_nu(ipp, im, ido);
// actual values of strains
long ncomp=Mm->ip[ipp].ncompstr;
strastrestate ssst = Mm->ip[ipp].ssst;
//vector epsp(ASTCKVEC(ncomp));
vector epst(ASTCKVEC(6));
vector eps(ASTCKVEC(ncomp));
vector epsp(ASTCKVEC(ncomp));
vector sig(ASTCKVEC(ncomp));
vector q(ASTCKVEC(1));
// actual total volumetric strain
Mm->givestrain(0, ipp, eps);
give_full_vector(epst, eps, ssst);
double epsv = first_invar(epst);
// actual plastic volumteric strain
//Mm->giveother(ipp, ido, ncomp, epsp);
// total volumetric strain from the last converged step
double epsvo = Mm->ip[ipp].eqother[ido+ncomp+7];
// total volumetric plastic strain from the last converged step
double epsvpo = Mm->ip[ipp].eqother[ido+ncomp+8];
// actual plastic strain
for(long i=0; i<ncomp; i++)
epsp(i) = Mm->ip[ipp].other[ido];
give_full_vector(epst, epsp, ssst);
double epsvp = first_invar(epst);
// actual value of p_c
//updateq(ipp, ido, eps, epsp, sig, q);
//double pc = q(0);
double pc = Mm->ip[ipp].other[ido+ncomp+1];
/*
// initial value of elastic strain
for (long i=0; i<ncomp; i++)
eps(i) = Mm->ic[ipp][3+i];
give_full_vector(auxt, eps, ssst);
double epsv_ini = first_invar(auxt);
if (epsv >= epsv_ini)
return e_ini;*/
// actual value of of specific volume
double v_ini = Mm->ip[ipp].eqother[ido+ncomp+4];
//double v = v_ini*(1.0+epsv);
//double p = pc/exp((v-v_pc)/kappa);
double po = Mm->ip[ipp].eqother[ido+ncomp+5];
if ((po == 0.0) || (epsv == 0.0))
return e_ini;
if (po > pc){
// for the unloading: \Delta eps^{el}_V = \Delta eps_V
// actual p according actual swelling line
double epsve = (epsv-epsvo) - (epsvp - epsvpo);
double p = po*exp(-(epsve)*v_ini/kappa);
if (epsve == 0.0)
e = 3.0*(1.0-2.0*nu)*v_ini*(-p)/kappa;
else
e = 3.0*(1.0-2.0*nu)*(p-po)/(epsve);
}
else{
e = e_ini;
}
//e = 3.0*(1.0-2.0*nu)*v*(-p)/kappa;
return e;
}
/**
The function returns actual Poisson's ratio.
@param ipp - integration point number
@param im - index of material type for given ip
@param ido - index of internal variables for given material in the ipp other array
*/
double dsmmat::give_actual_nu(long ipp, long im, long ido)
{
long idem = Mm->ip[ipp].gemid();
long ncompo = Mm->givencompeqother(ipp, im);
double nu = Mm->give_initial_nu(ipp, idem, ido+ncompo);
return nu;
}
/**
The function returns time step scaling factor required by the model. It is represented by ratio
of required time step size to the actual one or 1.0 for no change in time step size.
@param ipp - integration point number in the mechmat ip array.
@param im - index of material type for given ip
@param ido - index of internal variables for given material in the ipp other array
@return The function returns required time step scaling factor provided by the dsmmat model
Created by Tomas Koudelka, 5.2022
*/
double dsmmat::dstep_red(long ipp, long im, long ido)
{
long ncomp = Mm->ip[ipp].ncompstr;
return (Mm->ip[ipp].other[ido+ncomp+14]);
}
/**
Function returns irreversible plastic strains.
@param ipp - integration point number in the mechmat ip array.
@param ido - index of the first internal variable for given material in the ipp other array
@param epsp - %vector of irreversible strains
Returns vector of irreversible strains via parameter epsp
Created by Tomas Koudelka, 09.2012
*/
void dsmmat::giveirrstrains (long ipp, long ido, vector &epsp)
{
long i;
for (i=0;i<epsp.n;i++)
epsp[i] = Mm->ip[ipp].eqother[ido+i];
}
/**
This function extracts consistency parametr gamma for the reached equilibrium state
from the integration point other array.
@param ipp - integration point number in the mechmat ip array.
@param ido - index of internal variables for given material in the ipp other array
@retval The function returns value of consistency parameter.
Created by Tomas Koudelka, 09.2012
*/
double dsmmat::give_consparam (long ipp,long ido)
{
long ncompstr;
double gamma;
ncompstr=Mm->ip[ipp].ncompstr;
gamma = Mm->ip[ipp].eqother[ido+ncompstr];
return gamma;
}
/**
The function extracts saturation degree s from the integration point eqother array.
This value remains constant from initialization of the model.
@param ipp - integration point number in the mechmat ip array.
@param ido - index of internal variables for given material in the ipp other array
@retval The function returns value of saturation degree.
Created by Tomas Krejci, 14/12/2014
*/
double dsmmat::give_saturation_degree(long ipp, long ido)
{
long ncompstr;
double s;
ncompstr=Mm->ip[ipp].ncompstr;
s = Mm->ip[ipp].other[ido+ncompstr+9];
return s;
}
/**
The function extracts preconsolidation pressure p_c for the attained state
from the integration point other array.
@param ipp - integration point number in the mechmat ip array.
@param ido - index of internal variables for given material in the ipp other array
@retval The function returns value of preconsolidation pressure.
Created by Tomas Koudelka, 8.10.2013
*/
double dsmmat::give_preconspress(long ipp, long ido)
{
long ncompstr;
double pc;
ncompstr=Mm->ip[ipp].ncompstr;
pc = Mm->ip[ipp].other[ido+ncompstr+1];
return pc;
}
/**
The function extracts virgin void ratio e_lambda1 for the attained state (it is constant initial value)
from the integration point eqother array. This value remains constant from initialization of the model.
@param ipp - integration point number in the mechmat ip array.
@param ido - index of internal variables for given material in the ipp other array
@retval The function returns value of virgin porosity e_lambda1.
Created by Tomas Koudelka, 8.10.2013
*/
double dsmmat::give_virgporosity(long ipp, long ido)
{
long ncompstr;
double e_lambda1;
ncompstr=Mm->ip[ipp].ncompstr;
e_lambda1 = Mm->ip[ipp].eqother[ido+ncompstr+3]-1.0;
return e_lambda1;
}
/**
The function extracts initial porosity n_ini for the attained state (it is constant initial value)
from the integration point eqother array. This value remains constant from initialization of the model.
@param ipp - integration point number in the mechmat ip array.
@param ido - index of internal variables for given material in the ipp other array
@retval The function returns value of initial porosity e_ini.
Created by Tomas Koudelka, 8.10.2013
*/
double dsmmat::give_iniporosity(long ipp, long ido)
{
long ncompstr;
double v_ini, n_ini;
ncompstr=Mm->ip[ipp].ncompstr;
v_ini = Mm->ip[ipp].eqother[ido+ncompstr+4]; // v_ini = 1.0 + e_ini
n_ini = (v_ini-1.0)/v_ini;
return n_ini;
}
/**
The function extracts porosity for the actual state from the integration point other array.
@param ipp - integration point number in the mechmat ip array.
@param ido - index of internal variables for given material in the ipp other array
@retval The function returns value of actual porosity e.
Created by Tomas Koudelka, 06.2020
*/
double dsmmat::give_porosity(long ipp, long ido)
{
long ncompstr;
double n;
ncompstr=Mm->ip[ipp].ncompstr;
n = Mm->ip[ipp].other[ido+ncompstr+11];
return (n);
}
/**
The function returns rate of the volumetric strain at the given integration point.
@param ipp - integration point number in the mechmat ip array.
@param ido - index of internal variables for given material in the ipp other array
@return The function returns rate of the volumetric strain.
Created by Tomas Koudelka 06.2020
*/
double dsmmat::give_strain_vol_rate(long ipp, long ido)
{
long ncomp = Mm->ip[ipp].ncompstr;
return Mm->ip[ipp].other[ido+ncomp+10];
}
/**
The funtion marks required non-mechanical quantities in the array anmq.
@param anmq - array with flags for used material types
anmq[i] = 1 => qunatity type nonmechquant(i+1) is required
anmq[i] = 0 => qunatity type nonmechquant(i+1) is not required
@return The function does not return anything, but it may change content of anmq array.
*/
void dsmmat::give_reqnmq(long *anmq)
{
//if (pr_suc_f == no){
anmq[saturation_degree-1] = 1;
anmq[suction-1] = 1;
//}
//if (pr_tempr_f == no)
//anmq[temperature-1] = 1;
}
/**
The function changes material parameters for stochastic analysis.
@param atm - selected material parameters (parameters which are changed)
@param val - array containing new values of parameters
@return The function does not return anything.
Created by Tomas Koudelka, 09.2012
*/
void dsmmat::changeparam (atsel &atm,vector &val)
{
long i;
for (i=0;i<atm.num;i++)
{
switch (atm.atrib[i])
{
case 0:
m=val[i];
break;
case 1:
lambda=val[i];
break;
case 2:
kappa=val[i];
break;
case 3:
c1_tilde=val[i];
break;
case 4:
c2=val[i];
break;
case 5:
ks=val[i];
break;
default:
print_err("wrong number of atribute", __FILE__, __LINE__, __func__);
}
}
}

/trunk/SIFEL/MEFEL/SRC/dsm.h
0,0 → 1,148
#ifndef DSMMAT_H
#define DSMMAT_H
#include "galias.h"
#include "gfunct.h"
#include "xfile.h"
#include "alias.h"
#include "strretalg.h"
struct matrix;
struct vector;
struct atsel;
class nonlinman;
/**
This class defines double structure plasticity model for expansive clays based on cam-clay material model coupled with moisture transfer
and on double structure concept according to:
Sanchez, M., Gens, A., Do Nascimento Guimares, L., and Olivella, S., A double structure generalized plasticity model for expansive materials,
International Journal for Numerical and Analytical Methods in Geomechanics, vol. 29, no. 8, pp. 751-787, 2005. doi:10.1002/nag.434.
The material is defined by
three usual camclay material constants :
m - is frictional constant and it depends on the frictional angle phi
lambda - slope of the normal consolidationa line
kappa - slope of the swelling line
Additionally, there are three parameters which controls
evolution of hardening with respect to degree of saturation and suction
c1_tilde, c2 - controls the ratio of specific volumes v/v_sat (partially saturated/fully saturated)
ks - controls apparent adhesion (intersection of yield surface with p axis)
In addition to these material constants, several initial values
have to be specified at each integration point in the following order :
v_kappa1 - initial specific volume at the reference pressure p1 after
unloading from initial consolidation pressure p_c0 (v_kappa0)
p1 - reference pressure (p_1), it must be negative
p_c0 - initial consolidation pressure, it must be negative
sig0_1 -
. \
. \ components of eigenstresses
. /
sig0_ncomp -
Order of the eqother array components
id | description
------------------+------------------------------------------------
<0; ncompstr-1> | plastic strains
ncompstr | gamma - consistency parameter
ncompstr+1 | hardening parameter p_c - consolidation pressure
ncompstr+2 | hardening parameter p_s - apparent adhesion
ncompstr+3 | v_lambda1 - initial value of specific volume on the NCL at reference pressure p1
. | v_ini - initial specific volume
. | i1s - mean stress
. | j2s - the second invariant of stress deviator
. | epsv - total volume strain
ncompstr+8 | epsvp - plastic volume strain
ncompstr+9 | degree of saturation
ncompstr+10 | depsv/dt - volumtric strain rate
ncompstr+11 | actual porosity e
10.2012
TKo
*/
class dsmmat
{
public:
dsmmat (void);
~dsmmat (void);
void read (XFILE *in);
void initval(long ipp, long im, long ido);
double cohesion(vector &qtr);
double yieldfunction (vector &sig, vector &q);
void deryieldfsigma (vector &sig, vector &q, vector &dfds);
void dderyieldfsigma (matrix &ddfds);
void derpotsigma (vector &sig, vector &q, vector &dgds);
void deryieldfq(vector &sig, vector &q, vector &dq);
void deryieldfdqdq(matrix &ddfdq);
void deryieldfdsdq(matrix &dfdsdqt);
void der_q_gamma(long ipp, long ido, vector &sig, vector &qtr, vector &epsp, vector &dqdg);
void dqpardgammadsigma(long ipp, long ido, vector &sig, vector &qtr, vector &epsp, matrix &dqdgds);
void dqpardgammadqpar(long ipp, long ido, vector &sig, vector &qtr, vector &epsp, matrix &dqdgdq);
double plasmodscalar(long ipp, long ido, vector &sig, vector &qtr, vector &epsp);
void updateq(long ipp, long ido, vector &eps, vector &epsp, vector &sig, vector &q);
// void updateq(vector &epsp, strastrestate ssst, vector &q);
void matstiff (matrix &d, long ipp,long ido);
void depmatstiff (matrix &d, long ipp, long ido);
void nlstressesincr (long ipp,long im,long ido);
void givestressincr (long lcid, long ipp, long im, long ido, long fi, vector &sig);
void nlstresses (long ipp,long im,long ido);
long actualize_stiffmat(const nonlinman &nlman, long istep, long jstep);
void updateval (long ipp, long im, long ido);
double comp_actual_ym(long ipp, long im, long ido);
double give_actual_ym(long ipp, long im, long ido);
double give_actual_nu(long ipp, long im, long ido);
double dstep_red(long ipp, long im, long ido);
void giveirrstrains (long ipp, long ido, vector &epsp);
double give_consparam (long ipp, long ido);
double give_preconspress(long ipp, long ido);
double give_saturation_degree(long ipp, long ido);
double give_virgporosity(long ipp, long ido);
double give_iniporosity(long ipp, long ido);
double give_porosity(long ipp, long ido);
double give_strain_vol_rate(long ipp, long ido);
void give_reqnmq(long *anmq);
void changeparam (atsel &atm, vector &val);
/// frictional constant
double m;
/// slope of normal consolidation line
double lambda;
/// slope of swelling line
double kappa;
/// first parameter for control of the v/v_sat ratio evolution
double c1_tilde;
/// second parameter for control of the v/v_sat ratio evolution
double c2;
/// coefficent of apparent adhesion due to suction pressure
double ks;
/// soil swelling index reference value
double kappa_s0;
/// parameters for the strain rate caculation caused by the suction change
double aks1,aks2;
/// parameters for the thermal strain rate caculation
double a0,a2;
/// flag for suction definition (yes=suction evolution is prescribed by suc_fn, no=suction is taken from TRFEL)
//answertype pr_suc_f;
/// function of prescribed suction evolution
//gfunct suc_fn;
/// flag for temperature definition (yes=temperature is prescribed by tempr_fn, no=temperature is taken from TRFEL)
//answertype pr_tempr_f;
/// function of prescribed temperature evolution
//gfunct tempr_fn;
/// stress return algorithm
strretalg sra;
};
#endif

/trunk/SIFEL/MEFEL/SRC/elemswitch.cpp
5637,23 → 5637,25
#ifdef INC_OPENMP
#pragma omp parallel num_threads(Numth)
#pragma omp for
{
#pragma omp for
#endif
// just call computenlstresses procedure on all auxiliary integration points
// and check for math errors
for (app=0; app<Mm->tnaip; app++)
{
eid = Mm->elip[app];
if (Gtm->leso[eid] == 1)
{
Mm->computenlstresses(app,Mm->ip[app]);
err = check_math_errel(Mm->elip[app]);
if (err)
abort();
// just call computenlstresses procedure on all auxiliary integration points
// and check for math errors
for (app=0; app<Mm->tnaip; app++){
eid = Mm->elip[app];
if (Gtm->leso[eid] == 1){
Mm->computenlstresses(app,Mm->ip[app]);
err = check_math_errel(Mm->elip[app]);
if (err)
abort();
}
}
#ifdef INC_OPENMP
}
#endif
// swap regular integration point and auxiliary integration point arrays
// from now, all functions will work with REGULAR integration points
Mm->tnip = Mm->tnrip;

/trunk/SIFEL/MEFEL/SRC/globmat.cpp
60,7 → 60,6
abort ();
//FILE *ma;
//ma=fopen ("matice.txt","w");
513,8 → 512,8
// function computes strains at integration points
if (Mp->strcomp)
compute_ipstrains (lcid);
#ifdef DEBUGELEMFRC
#ifdef DEBUGELEMFRC
if (Elemfrc == NULL){
Elemfrc = new vector[Mt->ne];
for (i=0; i<Mt->ne; i++){
522,12 → 521,12
reallocv(ndofe, Elemfrc[i]);
}
}
#endif
#endif
// number of elements
long ne=Mt->ne;
//ivector seq(ne);
//long g = 0;
long ne=Mt->ne;
ivector seq(ne);
long g = 0;
#ifdef INC_OPENMP
double druntime = omp_get_wtime();
#else
538,9 → 537,9
#ifdef INC_OPENMP
double mstresscomp[6] = {}; // reduction array for integrated macrostress contributions from elements
#pragma omp parallel \
num_threads(Numth),\
private(i, j, ndofe, ndofemn, eid, cn, err, ifor, aifor),\
firstprivate(ncomp)
num_threads(Numth),\
private(i, j, ndofe, ndofemn, eid, cn, err, ifor, aifor),\
firstprivate(ncomp)
{
#pragma omp for reduction(+:mstresscomp[:6])
581,7 → 580,7
// this case means homogenization with some prescribed macro-stress components
if (ncomp == 0L)
ncomp = ndofemn-ndofe;
// localize contribution from standard element stresses due
// to macro-strains and fluctuations
// (macro-strains are added in compute_ipstrains function)
607,7 → 606,6
}
else{
// this case means hanging nodes
// the element contains hanging node
// the stiffness matrix has to be transformed
reallocv (RSTCKVEC(ndofemn,aifor));
617,18 → 615,18
}
else{
// the element does not contain hanging node nor homogenization
// #ifdef INC_OPENMP
// #pragma omp critical
// #endif
// {
#ifdef INC_OPENMP
#pragma omp critical
#endif
{
locglob (intfor,ifor.a,cn.a,ndofe);
//seq(i) = g++;
// }
seq(i) = g++;
}
}
}
#ifdef DEBUGELEMFRC
#ifdef DEBUGELEMFRC
copyv(ifor, Elemfrc[i]);
#endif
#endif
}
#ifdef INC_OPENMP
}
647,32 → 645,37
druntime = (clock() - druntime)/(double)CLOCKS_PER_SEC;
Omp_wtime += druntime;
#endif
// printv(seq, Out);
#ifdef DEBUGELEMFRC
fprintf(Out, "\nistep=%ld, jstep=%ld\n", Mp->istep, Mp->jstep);
// nullv (intfor,Ndofm);
#ifdef INC_OPENMP
//fprintf(Out, "Element sequence for istep=%ld, jstep=%ld:\n", Mp->istep, Mp->jstep);
//printv(seq, Out);
#endif
#ifdef DEBUGELEMFRC
//fprintf(Out, "\nistep=%ld, jstep=%ld\n", Mp->istep, Mp->jstep);
nullv (intfor,Ndofm);
for(i=0; i<Mt->ne; i++){
/* ndofe=Mt->give_ndofe(i);
ndofe=Mt->give_ndofe(i);
ndofemn = Gtm->give_ndofe (i);
reallocv (RSTCKIVEC(ndofemn,cn));
Mt->give_code_numbers (i,cn.a);
locglob (intfor,Elemfrc[i].a, cn.a, ndofe);*/
fprintf(Out, "eid=%4ld\n", i+1);
for(j=0; j<Elemfrc[i].n; j++)
fprintf(Out, " % .15le", Elemfrc[i][j]);
fprintf(Out, "\n");
locglob (intfor,Elemfrc[i].a, cn.a, ndofe);
//fprintf(Out, "eid=%4ld\n", i+1);
//for(j=0; j<Elemfrc[i].n; j++)
//fprintf(Out, " % .15le", Elemfrc[i][j]);
//fprintf(Out, "\n");
}
fprintf(Out, "\nElement order: istep=%ld, jstep=%ld\n", Mp->istep, Mp->jstep);
for(i=0; i<Mt->ne; i++)
fprintf(Out, "%5ld\n", h(i));
/*
fprintf(Out, "\nInternal force vector: istep=%ld, jstep=%ld\n", Mp->istep, Mp->jstep);
for(i=0; i<Ndofm; i++){
fprintf(Out, " % .15le\n", intfor[i]);
}*/
fflush(Out);
#endif
//fprintf(Out, "\nElement order: istep=%ld, jstep=%ld\n", Mp->istep, Mp->jstep);
//for(i=0; i<Mt->ne; i++)
//fprintf(Out, "%5ld\n", h(i));
//fprintf(Out, "\nInternal force vector: istep=%ld, jstep=%ld\n", Mp->istep, Mp->jstep);
//for(i=0; i<Ndofm; i++){
//fprintf(Out, " % .15le\n", intfor[i]);
//}
//fflush(Out);
#endif
}

/trunk/SIFEL/MEFEL/SRC/homogmatm.cpp
0,0 → 1,414
/*
File: homogmatm.cpp
Author: Tomas Krejci, 22/07/2022
Purpose: Calculates properties of homogenized elastic material
*/
#include "mechmat.h"
#include "probdesc.h"
#include "homogmatm.h"
#include "matrix.h"
#include "vector.h"
#include "stochdriver.h"
#include "global.h"
#include "intpoints.h"
#include "vecttens.h"
/**
Constructor initializes data members to zero or default values.
Created by TKr,
*/
homogmatm::homogmatm (void)
{
allocm (6,6,dd);
hom_mattypem = 0;
hom_mattypem_number = 0;
}
/**
Destructor is defined only for the formal purposes.
Created by TKr,
*/
homogmatm::~homogmatm (void)
{
destrm(dd);
}
/**
Function reads material parameters from the opened text file.
@param in - pointer to the opened XFILE
@return The function does not return anything.
Created by TKr,
*/
void homogmatm::read (XFILE *in)
{
xfscanf (in,"%ld %ld",&hom_mattypem,&hom_mattypem_number);
hom_mattypem_number = hom_mattypem_number-1;
}
/**
Function prints material parameters into the opened text file.
@param out - pointer to the opened FILE
@return The function does not return anything.
Created by TKr
*/
void homogmatm::print (FILE *out)
{
fprintf (out," %ld %ld",hom_mattypem,hom_mattypem_number+1);
}
/**
Function assembles stiffness %matrix of material.
@param d - stiffness %matrix of material (output)
@param ssst - strain/stress state
@return The function returns material stiffness %matrix in the parameter d.
Created by TKr
*/
void homogmatm::matstiff (matrix &d,strastrestate ssst)
{
switch (ssst){
case planestress:{
matstiff_plstress (d);
break;
}
case planestrain:{
matstiff_plstrain (d);
break;
}
case axisymm:{
matstiff_axi (d);
break;
}
case spacestress:{
matstiff_spacestr (d);
break;
}
default:{
print_err("unknown number of components of stress tensor is required", __FILE__, __LINE__, __func__);
}
}
}
/**
Function creates stiffness matrix of the elastic
isotropic material for 2D problems (plane stress).
@param d - stiffness matrix of the material
@return The function returns material stiffness %matrix in the parameter d.
Created by TKr
*/
void homogmatm::matstiff_plstress (matrix &d)
{
d[0][0] = dd[0][0]; d[0][1] = dd[0][1]; d[0][2] = dd[0][2];
d[1][0] = dd[1][0]; d[1][1] = dd[1][1]; d[1][2] = dd[1][2];
d[2][0] = dd[2][0]; d[2][1] = dd[2][1]; d[2][2] = dd[2][2];
}
/**
Function creates stiffness %matrix of the elastic
isotropic material for 2D problems (plane strain).
@param d - stiffness %matrix of the material
@return The function returns material stiffness %matrix in the parameter d.
Created by TKr
*/
void homogmatm::matstiff_plstrain (matrix &d)
{
d[0][0] = dd[0][0]; d[0][1] = dd[0][1]; d[0][2] = dd[0][2];
d[1][0] = dd[1][0]; d[1][1] = dd[1][1]; d[1][2] = dd[1][2];
d[2][0] = dd[2][0]; d[2][1] = dd[2][1]; d[2][2] = dd[2][2];
// zde rozmyslet
if (d.m > 3)
{
d[0][3] = dd[0][1]; d[1][3] = dd[1][0];
d[3][0] = dd[1][0]; d[3][1] = dd[1][0]; d[3][3] = dd[1][1];
}
}
/**
Function creates stiffness %matrix of the elastic
isotropic material for 2D axisymmetric problems.
@param d - stiffness %matrix of the material
@return The function returns material stiffness %matrix in the parameter d.
Created by TKr, 19.7.2001
*/
void homogmatm::matstiff_axi (matrix &d)
{
d[0][0] = dd[0][0]; d[0][1] = dd[0][1]; d[0][2] = dd[0][2]; d[0][3] = dd[0][3];
d[1][0] = dd[1][0]; d[1][1] = dd[1][1]; d[1][2] = dd[1][2]; d[1][3] = dd[1][3];
d[2][0] = dd[2][0]; d[2][1] = dd[2][1]; d[2][2] = dd[2][2]; d[2][3] = dd[2][3];
d[3][0] = dd[3][0]; d[3][1] = dd[3][1]; d[3][2] = dd[3][2]; d[3][3] = dd[3][3];
}
/**
Function creates stiffness %matrix of the elastic
isotropic material for 3D problems.
@param d - stiffness %matrix of the material
Created by TKr
*/
void homogmatm::matstiff_spacestr (matrix &d)
{
d[0][0] = dd[0][0]; d[0][1] = dd[0][1]; d[0][2] = dd[0][2]; d[0][3] = dd[0][3]; d[0][4] = dd[0][4]; d[0][5] = dd[0][5];
d[1][0] = dd[1][0]; d[1][1] = dd[1][1]; d[1][2] = dd[1][2]; d[1][3] = dd[1][3]; d[1][4] = dd[1][4]; d[1][5] = dd[1][5];
d[2][0] = dd[2][0]; d[2][1] = dd[2][1]; d[2][2] = dd[2][2]; d[2][3] = dd[2][3]; d[2][4] = dd[2][4]; d[2][5] = dd[2][5];
d[3][0] = dd[3][0]; d[3][1] = dd[3][1]; d[3][2] = dd[3][2]; d[3][3] = dd[3][3]; d[3][4] = dd[3][4]; d[3][5] = dd[3][5];
d[4][0] = dd[4][0]; d[4][1] = dd[4][1]; d[4][2] = dd[4][2]; d[4][3] = dd[4][3]; d[4][4] = dd[4][4]; d[4][5] = dd[4][5];
d[5][0] = dd[5][0]; d[5][1] = dd[5][1]; d[5][2] = dd[5][2]; d[5][3] = dd[5][3]; d[5][4] = dd[5][4]; d[5][5] = dd[5][5];
}
/**
Function assembles compliance %matrix of material.
@param c - complience %matrix of material (output)
@param ssst - strain/stress state
@return The function returns material complience %matrix in the parameter d.
Created by TKr
*/
void homogmatm::matcompl (matrix &/*c*/,strastrestate ssst)
{
matrix d;
allocm (6,6,d);
matstiff (d,ssst);
//c = invm(d);
destrm(d);
}
/**
function assembles the conductivity %matrix after homogenization
@param d - array of %matrix entries
Created by TKr
*/
void homogmatm::assemble_matrices (double *d,long ncomp,long /*dim*/)
{
if (ncomp == 3){
dd[0][0]=d[0]; dd[0][1]=d[1]; dd[0][2]=d[2];
dd[1][0]=d[3]; dd[1][1]=d[4]; dd[1][2]=d[5];
dd[2][0]=d[6]; dd[2][1]=d[7]; dd[2][2]=d[8];
/* fprintf (Out,"\n");
fprintf (Out,"\n\n d[0] = %e",d[0]);
fprintf (Out,"\n\n d[1] = %e",d[1]);
fprintf (Out,"\n\n d[2] = %e",d[2]);
fprintf (Out,"\n\n d[3] = %e",d[3]);
fprintf (Out,"\n\n d[4] = %e",d[4]);
fprintf (Out,"\n\n d[5] = %e",d[5]);
fprintf (Out,"\n\n d[6] = %e",d[6]);
fprintf (Out,"\n\n d[7] = %e",d[7]);
fprintf (Out,"\n\n d[8] = %e",d[8]);
fprintf (Out,"\n");
fflush(Out);
*/
}
if (ncomp == 4){
dd[0][0]=d[0]; dd[0][1]=d[1]; dd[0][2]=d[2]; dd[0][3]=d[3];
dd[1][0]=d[4]; dd[1][1]=d[5]; dd[1][2]=d[6]; dd[1][3]=d[7];
dd[2][0]=d[8]; dd[2][1]=d[9]; dd[2][2]=d[10]; dd[2][3]=d[11];
dd[3][0]=d[12]; dd[3][1]=d[13]; dd[3][2]=d[14]; dd[3][3]=d[15];
}
if (ncomp == 6){
dd[0][0]=d[0]; dd[0][1]=d[1]; dd[0][2]=d[2]; dd[0][3]=d[3]; dd[0][4]=d[4]; dd[0][5]=d[5];
dd[1][0]=d[6]; dd[1][1]=d[7]; dd[1][2]=d[8]; dd[1][3]=d[9]; dd[1][4]=d[10]; dd[1][5]=d[11];
dd[2][0]=d[12]; dd[2][1]=d[13]; dd[2][2]=d[14]; dd[2][3]=d[15]; dd[2][4]=d[16]; dd[2][5]=d[17];
dd[3][0]=d[18]; dd[3][1]=d[19]; dd[3][2]=d[20]; dd[3][3]=d[21]; dd[3][4]=d[22]; dd[3][5]=d[23];
dd[4][0]=d[24]; dd[4][1]=d[25]; dd[4][2]=d[26]; dd[4][3]=d[27]; dd[4][4]=d[28]; dd[4][5]=d[29];
dd[5][0]=d[30]; dd[5][1]=d[31]; dd[5][2]=d[32]; dd[5][3]=d[33]; dd[5][4]=d[34]; dd[5][5]=d[35];
/* fprintf (Out,"\n");
fprintf (Out,"\n\n d[0] = %e",d[0]);
fprintf (Out,"\n\n d[1] = %e",d[1]);
fprintf (Out,"\n\n d[2] = %e",d[2]);
fprintf (Out,"\n\n d[3] = %e",d[3]);
fprintf (Out,"\n\n d[4] = %e",d[4]);
fprintf (Out,"\n\n d[5] = %e",d[5]);
fprintf (Out,"\n\n d[6] = %e",d[6]);
fprintf (Out,"\n\n d[7] = %e",d[7]);
fprintf (Out,"\n\n d[8] = %e",d[8]);
fprintf (Out,"\n\n d[9] = %e",d[9]);
fprintf (Out,"\n\n d[10] = %e",d[10]);
fprintf (Out,"\n\n d[11] = %e",d[11]);
fprintf (Out,"\n\n d[12] = %e",d[12]);
fprintf (Out,"\n\n d[13] = %e",d[13]);
fprintf (Out,"\n\n d[14] = %e",d[14]);
fprintf (Out,"\n\n d[15] = %e",d[15]);
fprintf (Out,"\n\n d[16] = %e",d[16]);
fprintf (Out,"\n\n d[17] = %e",d[17]);
fprintf (Out,"\n\n d[18] = %e",d[18]);
fprintf (Out,"\n\n d[19] = %e",d[19]);
fprintf (Out,"\n\n d[20] = %e",d[20]);
fprintf (Out,"\n\n d[21] = %e",d[21]);
fprintf (Out,"\n\n d[22] = %e",d[22]);
fprintf (Out,"\n\n d[23] = %e",d[23]);
fprintf (Out,"\n\n d[24] = %e",d[24]);
fprintf (Out,"\n\n d[25] = %e",d[25]);
fprintf (Out,"\n\n d[26] = %e",d[26]);
fprintf (Out,"\n\n d[27] = %e",d[27]);
fprintf (Out,"\n\n d[28] = %e",d[28]);
fprintf (Out,"\n\n d[29] = %e",d[29]);
fprintf (Out,"\n\n d[30] = %e",d[30]);
fprintf (Out,"\n\n d[31] = %e",d[31]);
fprintf (Out,"\n\n d[32] = %e",d[32]);
fprintf (Out,"\n\n d[33] = %e",d[33]);
fprintf (Out,"\n\n d[34] = %e",d[34]);
fprintf (Out,"\n\n d[35] = %e",d[35]);
fprintf (Out,"\n");
fflush(Out);
*/
}
}
/**
This function initializes material model data
with respect of consistency parameter gamma.
@param ipp - integration point number
@param ido - index of internal variables for given material in the ipp other array
Created by TKr
*/
void homogmatm::initval(long /*ipp*/, long /*im*/, long /*ido*/)
{
//Mm->ip[ipp].other[ido+0] = Mm->ip[ipp].eqother[ido+0] = 0.0;
//Mm->ip[ipp].other[ido+1] = Mm->ip[ipp].eqother[ido+1] = 0.0;
}
/**
Function computes true stresses.
@param ipp - number of integration point
@return The function does not return anything.
Created by TKr
*/
void homogmatm::nlstresses (long ipp, long /*im*/, long /*ido*/)
{
long i, n = Mm->ip[ipp].ncompstr;
vector eps(n),sig(n),epst(ASTCKVEC(6));
strastrestate ssst = Mm->ip[ipp].ssst;
matrix d(n,n);
// initial values
for (i=0;i<n;i++){
eps[i]=Mm->ip[ipp].strain[i];
}
matstiff (d,ssst);
mxv (d,eps,sig);
/* //it is necessary to complete here:
if (Mm->ip[ipp].ssst == planestress)
Mm->ip[ipp].strain[3] = -nu / (1.0 - nu) * (eps[0]+eps[1]);
*/
if ((Mp->eigstrains == 4) || (Mp->eigstrains == 5))
{
for (i=0;i<n;i++)
sig(i) += Mm->eigstresses[ipp][i];
}
for (i=0;i<n;i++){
Mm->ip[ipp].stress[i]=sig(i);
}
}
/**
This function updates values in the other array reached in the previous equlibrium state to
values reached in the new actual equilibrium state.
@param ipp - integration point number in the mechmat ip array.
@param ido - index of internal variables for given material in the ipp other array
Created by TKr
*/
void homogmatm::updateval (long /*ipp*/,long /*ido*/)
{
//Mm->ip[ipp].eqother[ido+0] = Mm->ip[ipp].other[ido+0]; //
//Mm->ip[ipp].eqother[ido+1] = Mm->ip[ipp].other[ido+1]; //
}
/**
The function returns rate of the volumetric strain rate at the given integration point.
@param ipp - integration point number in the mechmat ip array.
@param ido - index of internal variables for given material in the ipp other array
@return The function returns rate of the volumetric strain.
Created by TKr according to Tomas Koudelka
*/
double homogmatm::give_strain_vol(long /*ipp*/, long /*ido*/)
{
//return Mm->ip[ipp].eqother[ido+1];
print_err("the function was not defined properly!", __FILE__, __LINE__, __func__);
return 0.0;
}

/trunk/SIFEL/MEFEL/SRC/homogmatm.h
0,0 → 1,49
#ifndef HOMOGMATM_H
#define HOMOGMATM_H
#include "alias.h"
#include "genfile.h"
#include "iotools.h"
struct matrix;
struct vector;
struct atsel;
/**
The class homogmatm defines elastic isotropic material model.
Material parameters are:
- Young's modulus (modulus of elasticity)
- Poisson's ratio
Created by TKr,
*/
class homogmatm
{
public:
homogmatm (void);
~homogmatm (void);
void read (XFILE *in);
void print (FILE *out);
void matstiff (matrix &d,strastrestate ssst);
void matstiff_plstress (matrix &d);
void matstiff_plstrain (matrix &d);
void matstiff_axi (matrix &d);
void matstiff_spacestr (matrix &d);
void matcompl (matrix &c,strastrestate ssst);
void assemble_matrices (double *d,long ncomp, long dim);
void initval(long ipp, long im, long ido);
void nlstresses (long ipp, long im, long ido);
void updateval (long ipp,long ido);
double give_strain_vol(long ipp, long ido);
matrix dd;
long hom_mattypem; //material type for mogenization
long hom_mattypem_number; //number of material for homogenization
};
#endif

/trunk/SIFEL/MEFEL/SRC/plquadiface.cpp
391,8 → 391,8
transf[1][0]=-dy/l;
transf[1][1]=dx/l;
// nodes 1 and 4 in local element ordering has natural coordinate ksi=-1
// nodes 2 and 3 in local element ordering has natural coordinate ksi=+1
// nodes 1 and 4 in local element ordering have natural coordinate ksi=-1
// nodes 2 and 3 in local element ordering have natural coordinate ksi=+1
// node 1 coincides with node 4
// node 2 coincides with node 3