WebSVN - sifel - Path Comparison - / Rev 775 and / Rev 776

Ignore whitespace Rev 775 → Rev 776

 /trunk/SIFEL/GEFEL/elemtools.cpp
 ,5 → 1,6
 #include <stdlib.h>
 #include <math.h>
+#include "genfile.h"
 #include "elemtools.h"
 #include "matrix.h"
 ,3 → 366,294
   return area;
 }
+/**
+   function computes volume of an hexahedral finite element
+   @param eid - element id (it is used only in the case of nonpositive volume)
+   @param x,y,z - arrays of node coordinates
+   TKr 26/09/2022 according to JK
+*/
+double hexahedra_volume (long eid,vector &x,vector &y,vector &z)
+{
+  long i,j,k;
+  double xi,eta,zeta,ww1,ww2,ww3,jac,vol;
+  vector w(2),gp(2);
+  gauss_points (gp.a,w.a,2);
+  vol=0.0;
+  for (i=0;i<2;i++){
+    xi=gp[i];  ww1=w[i];
+    for (j=0;j<2;j++){
+      eta=gp[j];  ww2=w[j];
+      for (k=0;k<2;k++){
+        zeta=gp[k]; ww3=w[k];
+        jac_3d (jac,x,y,z,xi,eta,zeta);
+        if (jac<0.0){
+          print_err("\n negative Jacobian on element %ld", __FILE__, __LINE__, __func__,eid);
+          jac=fabs(jac);
+        }
+        jac*=ww1*ww2*ww3;
+        vol+=jac;
+      }
+    }
+  }
+  return vol;
+}
+/**
+   function generates outer normal %vector to required element surface
+   @param surfid - id of surface
+   @param x,y,z - arrays of node coordinates
+   @param n - %vector of the outer normal
+   TKr 26/09/2022 according to JK
+*/
+void hexahedra_normal_vectors (long surfid,vector &x,vector &y,vector &z,vector &n)
+{
+  double norm;
+  double x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4;
+  vector t1(3),t2(3);
+  //two vectors are selected - one edge and one diagonal
+  if (surfid==0){
+    x1 = x[0];
+    x2 = x[3];
+    x3 = x[7];
+    x4 = x[4];
+    y1 = y[0];
+    y2 = y[3];
+    y3 = y[7];
+    y4 = y[4];
+    z1 = z[0];
+    z2 = z[3];
+    z3 = z[7];
+    z4 = z[4];
+  }
+  if (surfid==1){
+    x1 = x[1];
+    x2 = x[0];
+    x3 = x[4];
+    x4 = x[5];
+    y1 = y[1];
+    y2 = y[0];
+    y3 = y[4];
+    y4 = y[5];
+    z1 = z[1];
+    z2 = z[0];
+    z3 = z[4];
+    z4 = z[5];
+  }
+  if (surfid==2){
+    x1 = x[2];
+    x2 = x[1];
+    x3 = x[5];
+    x4 = x[6];
+    y1 = y[2];
+    y2 = y[1];
+    y3 = y[5];
+    y4 = y[6];
+    z1 = z[2];
+    z2 = z[1];
+    z3 = z[5];
+    z4 = z[6];
+  }
+  if (surfid==3){
+    x1 = x[3];
+    x2 = x[2];
+    x3 = x[6];
+    x4 = x[7];
+    y1 = y[3];
+    y2 = y[2];
+    y3 = y[6];
+    y4 = y[7];
+    z1 = z[3];
+    z2 = z[2];
+    z3 = z[6];
+    z4 = z[7];
+  }
+  if (surfid==4){
+    x1 = x[0];
+    x2 = x[1];
+    x3 = x[2];
+    x4 = x[3];
+    y1 = y[0];
+    y2 = y[1];
+    y3 = y[2];
+    y4 = y[3];
+    z1 = z[0];
+    z2 = z[1];
+    z3 = z[2];
+    z4 = z[3];
+  }
+  if (surfid==5){
+    x1 = x[4];
+    x2 = x[5];
+    x3 = x[6];
+    x4 = x[7];
+    y1 = y[4];
+    y2 = y[5];
+    y3 = y[6];
+    y4 = y[7];
+    z1 = z[4];
+    z2 = z[5];
+    z3 = z[6];
+    z4 = z[7];
+  }
+  // first in-plane vector
+  t1[0]=x2-x1;
+  t1[1]=y2-y1;
+  t1[2]=z2-z1;
+  //  second in-plane vector
+  t2[0]=x3-x1;
+  t2[1]=y3-y1;
+  t2[2]=z3-z1;
+  //  outer normal vector to the first element surface
+  n[0]=t1[1]*t2[2]-t1[2]*t2[1];
+  n[1]=t1[2]*t2[0]-t1[0]*t2[2];
+  n[2]=t1[0]*t2[1]-t1[1]*t2[0];
+  //  length of the normal vector
+  norm = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
+  if (norm<1.0e-10){
+    print_err("nonpositive length of the normal vector",__FILE__,__LINE__,__func__);
+    abort ();
+  }
+  n[0]/=norm;
+  n[1]/=norm;
+  n[2]/=norm;
+}
+/**
+   function computes areas of the tetrahedron surfaces
+   @param surfid - id of surface
+   @param x,y,z - arrays of node coordinates
+   TKr 26/09/2022 according to JK
+*/
+double hexahedra_surface_areas (long surfid,vector &x,vector &y,vector &z)
+{
+  double area,l,l1,l2,l3,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4;
+  if (surfid==0){
+    x1 = x[0];
+    x2 = x[3];
+    x3 = x[7];
+    x4 = x[4];
+    y1 = y[0];
+    y2 = y[3];
+    y3 = y[7];
+    y4 = y[4];
+    z1 = z[0];
+    z2 = z[3];
+    z3 = z[7];
+    z4 = z[4];
+  }
+  if (surfid==1){
+    x1 = x[1];
+    x2 = x[0];
+    x3 = x[4];
+    x4 = x[5];
+    y1 = y[1];
+    y2 = y[0];
+    y3 = y[4];
+    y4 = y[5];
+    z1 = z[1];
+    z2 = z[0];
+    z3 = z[4];
+    z4 = z[5];
+  }
+  if (surfid==2){
+    x1 = x[2];
+    x2 = x[1];
+    x3 = x[5];
+    x4 = x[6];
+    y1 = y[2];
+    y2 = y[1];
+    y3 = y[5];
+    y4 = y[6];
+    z1 = z[2];
+    z2 = z[1];
+    z3 = z[5];
+    z4 = z[6];
+  }
+  if (surfid==3){
+    x1 = x[3];
+    x2 = x[2];
+    x3 = x[6];
+    x4 = x[7];
+    y1 = y[3];
+    y2 = y[2];
+    y3 = y[6];
+    y4 = y[7];
+    z1 = z[3];
+    z2 = z[2];
+    z3 = z[6];
+    z4 = z[7];
+  }
+  if (surfid==4){
+    x1 = x[0];
+    x2 = x[1];
+    x3 = x[2];
+    x4 = x[3];
+    y1 = y[0];
+    y2 = y[1];
+    y3 = y[2];
+    y4 = y[3];
+    z1 = z[0];
+    z2 = z[1];
+    z3 = z[2];
+    z4 = z[3];
+  }
+  if (surfid==5){
+    x1 = x[4];
+    x2 = x[5];
+    x3 = x[6];
+    x4 = x[7];
+    y1 = y[4];
+    y2 = y[5];
+    y3 = y[6];
+    y4 = y[7];
+    z1 = z[4];
+    z2 = z[5];
+    z3 = z[6];
+    z4 = z[7];
+  }
+  //the quad (1,2,3,4) is split into two triangles
+  //the first triangle with nodes 1,2,4
+  l1=sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) + (z2-z1)*(z2-z1));
+  l2=sqrt((x4-x2)*(x4-x2) + (y4-y2)*(y4-y2) + (z4-z2)*(z4-z2));
+  l3=sqrt((x1-x4)*(x1-x4) + (y1-y4)*(y1-y4) + (z1-z4)*(z1-z4));
+  // Heron's formula for the first triangle area
+  l=(l1+l2+l3)/2.0;
+  area = sqrt(l*(l-l1)*(l-l2)*(l-l3));
+  //second triangle with nodes 2,3,4
+  l1=sqrt((x3-x2)*(x3-x2) + (y3-y2)*(y3-y2) + (z3-z2)*(z3-z2));
+  l2=sqrt((x4-x3)*(x4-x3) + (y4-y3)*(y4-y3) + (z4-z3)*(z4-z3));
+  l3=sqrt((x2-x4)*(x2-x4) + (y2-y4)*(y2-y4) + (z2-z4)*(z2-z4));
+  // Heron's formula for the triangle area
+  l=(l1+l2+l3)/2.0;
+  area = area + sqrt(l*(l-l1)*(l-l2)*(l-l3));
+  return area;
+}

 /trunk/SIFEL/GEFEL/elemtools.h
 ,6 → 27,12
 ///  function computes areas of the tetrahedron surfaces
 double tetrahedra_surface_areas (long surfid,vector &x,vector &y,vector &z);
+///  function computes volume of an element
+double hexahedra_volume (long eid,vector &x,vector &y,vector &z);
+///  function constructs the outer unit normal %vector to an element surface
+void hexahedra_normal_vectors (long surfid,vector &x,vector &y,vector &z,vector &n);
+///  function computes areas of the hexahedron surfaces
+double hexahedra_surface_areas (long surfid,vector &x,vector &y,vector &z);
 #endif

/trunk/SIFEL/GEFEL/gtopology.cpp
5906,6 → 5906,7

// reading of the array eldom
// this line is avoided only for purposes of Charles bridge
// zde rozmyslet??!!
stop->read_eldom (ne,in);
}

 /trunk/SIFEL/GEFEL/matrix.cpp
 ,16 → 280,18
 */
 long allocm(long m, long n, double **&mat)
 {
+  long i;
   mat = new double* [m];
   if (mat)
-  {
-    for (register long i=0; i<m; i++)
-      mat[i] = new double [n];
-    return (0);
-  }
+    {
+      for (i=0; i<m; i++)
+        mat[i] = new double [n];
+      return (0);
+    }
   else
     print_err("cannot allocate memory for array=matrix a(%ld,%ld)", __FILE__, __LINE__, __func__, m, n);
   return (1);
 }

 /trunk/SIFEL/HPARAL/SRC/Makefile
 ,19 → 1,20
-####################################################################
+################################################################################
 #
-#  Makefile for hybrid parallel mechanical and transport code HPARAL
+#  Makefile for hybrid parallel mechanical, transport, and coupled  code HPARAL
 #
-####################################################################
+################################################################################
 # Default values for make command line targets
 PARAM    = $(MAKECMDGOALS)
+TGTHPMET = -f Makefilec
 TGTHPMEF = -f Makefilem
 TGTHPTRF = -f Makefilet
 PARDEP = all
-ifneq ($(and $(findstring hpmefel,$(MAKECMDGOALS)), $(findstring hptrfel,$(MAKECMDGOALS))),)
+ifneq ($(and $(findstring hpmetr,$(MAKECMDGOALS)),  $(findstring hpmefel,$(MAKECMDGOALS)), $(findstring hptrfel,$(MAKECMDGOALS))),)
 #
-# set parameters in the case that hpmefel and hptrfel were both specified on the command line
+# set parameters in the case that hpmetr, hpmefel, and hptrfel were specified on the command line
 #
 PARAM = $(filter-out hpmefel,$(MAKECMDGOALS)) # remove 'hpmefel' from command line string
 ,7 → 21,8
 AUX = $(filter-out hptrfel,$(PARAM))          # remove 'hptrfel' from command line string
 PARAM := $(AUX)                               # store remaining target options
 PARDEP = ; @(echo)                            # avoiding make message 'Nothing to be done for...'
-TRFDEP = |hpmefel;                            # if hpmefel and hptrfel are specified, hptrfel waits on finish of hpmefel
+TRFDEP = |hpmefel;                            # if hpmefel and hptrfel are specified, hptrfel waits on finish of hpmefel and hpmetr waits on finish of hptrfel
+METDEP = |hptrfel;                            # if hpmefel and hptrfel are specified, hptrfel waits on finish of hpmefel and hpmetr waits on finish of hptrfel
 else
 #
 ,9 → 47,18
 endif
 endif
+# hpmetr was specified on the command line
+ifeq ($(findstring hpmetr,$(MAKECMDGOALS)), hpmetr)
+TGTHPMET =
+PARAM = $(filter-out hpmetr,$(MAKECMDGOALS)) # remove hpmetr from make parameters and left the remaining target options
+ifneq ($(PARAM),)
+PARDEP = ; @(echo)                            # avoiding the make message 'Nothing to be done for...'
 endif
+endif
+endif
 #########################################################################
 #  Targets
 #
 ,10 → 69,10
 .PHONY : cleanomp cleanompdeb cleanompopt cleanompoptdeb
 .PHONY : cleandep cleandepdeb clenadepopt cleandepoptdeb
 .PHONY : cleandepompdeb cleandepompopt cleandepompoptdeb
-.PHONY : hpmefel hptrfel
+.PHONY : hpmefel hptrfel hpmetr
-all: hpmefel hptrfel
+all: hpmefel hptrfel hpmetr
 hpmefel:
 ,7 → 88,13
 endif
+hpmetr: $(METDEP)
+ifneq ($(TGTHPMET),)
+        +@($(MAKE) $(TGTHPMET) $(PARAM))
+endif
 opt:    $(PARDEP)
 deb:    $(PARDEP)
 optdeb: $(PARDEP)

 /trunk/SIFEL/HPARAL/SRC/Makefilem
 ,7 → 36,7
 ####################################################################
 # List of source files
 #
-SRCS = hplssolver.cpp hpmefel.cpp hpmefelinit.cpp hpsolverm.cpp
+SRCS = hplssolver.cpp hpmefel.cpp hpmefelinit.cpp hpmtsolver.cpp hpsolver.cpp hpsolverm.cpp

 /trunk/SIFEL/HPARAL/SRC/hpglobal.h
 ,11 → 9,10
 #include "genfile.h"
 #include "seqfilesm.h"
 #include "seqfilest.h"
+#include "global.h"
 #include "globalt.h"
 //**************************************
 //**************************************
 //  GLOBAL CONSTANTS AND VARIABLES
 ,7 → 58,6
 EXTERN long** Tiling;
 //  type of problem on the MEFEL master processor
 //  it is needed also on slave processors
 EXTERN problemtype Mtprobm;

 /trunk/SIFEL/HPARAL/SRC/hplssolver.cpp
 ,7 → 14,7
    JK, 4. 6. 2013
 */
-void parallel_homogenization_linear_statics ()
+void parallel_homogenization_linear_statics_tiles ()
 {
   long i,j,combmin,combmax,ncomb,tid;
   double *lhs,*rhs,**condm,**condvlhs,**condvrhs;
 ,3 → 208,598
+/**
+   function computes linear static problem on a macro-level
+   material parameters are obtained by homogenization
+   on microscale, stationary problems are solved
+   these microproblems can be defined on every finite element or
+   on aggregates of finite elements, the aggregates can send
+   data from all aggregated elements or averaged data which
+   has the same form as data from a single element
+   material parameters are not defined in material models, they are
+   obtained from a meso or micro level by a homogenization method
+   TKr, 14/07/2022
+*/
+void parallel_homogenization_linear_statics ()
+{
+  long i,j,k,ii,jj,kk,n,aux,nelidom,buffsize1,buffsize2,ncomp,dim;
+  long *neldom,*buff;
+  double end_time,*lhs,*rhs,*fl,*fi;
+  double *buff1,*buff2;
+  MPI_Status stat;
+  //  allocation of buffer
+  //  buff[0] - the number of connections between the macroproblem and microproblem
+  //  buff[1] - the maximum number of connections between the macroproblem and microproblem
+  buff = new long [2];
+  if (Myrank==0){
+    switch (Mp->mami){
+    case elem_domain:{
+      //  each element is connected with its own microproblem
+      buff[0]=1;
+      buff[1]=1;
+      break;
+    }
+    case aggregate_domain:{
+      //  generation of auxiliary data needed in parallel computation
+      //  elements of macroproblem are aggregated, each finite element
+      //  sends its data to microproblem and obtains its matrices
+      //  map between finite elements of macroproblem and microproblems
+      //  this map is read in the subroutine transtop::read_seq_top
+      //  eldom[i]=j - the i-th finite element of macroproblem is connected to the j-th microproblem
+      //  the number of elements on subdomains
+      //  it is the number of macroelements connected to one microproblem
+      //  the number of subdomains--microproblems must be Nproc-1
+      //  neldom[i]=j - j finite elements of the macroproblem are connected to the i-th microproblem
+      //  neldom[0]=0 - the master processor (myrank=0) contains no microproblem
+      neldom = new long [Nproc];
+      for (i=0;i<Nproc;i++){
+        neldom[i]=0;
+      }
+      for (i=0;i<Mt->ne;i++){
+        neldom[Gtm->stop->eldom[i]]++;
+      }
+      buffsize1=0;
+      for (i=0;i<Nproc;i++){
+        if (buffsize1<neldom[i])
+          buffsize1=neldom[i];
+      }
+      buff[1]=buffsize1;
+      break;
+    }
+    case material_aggregate_domain:{
+      //one domain represents one material = one processor
+      //the aggregate covers elements of one materyal type
+      neldom = new long [Nproc];
+      for (i=0;i<Nproc;i++){
+        neldom[i]=0;
+      }
+      for (i=0;i<Mt->ne;i++)
+        neldom[Gtm->stop->eldom[i]] = 1;
+      buff[0]=1;
+      buff[1]=1;
+      break;
+    }
+    case eff_aggregate_domain:{
+      neldom = new long [Nproc];
+      for (i=0;i<Nproc;i++){
+        neldom[i]=0;
+      }
+      for (i=0;i<Mt->ne;i++){
+        if(Gtm->stop->eldom[i] != Mt->elements[i].idm[0])
+          neldom[Gtt->stop->eldom[i]]++;
+      }
+      //  the aggregate sends averaged data to microproblem
+      //  only aggregate sends and receives data, the same
+      //  matrices are used on all elements in the aggregate
+      buff[0]=1;
+      buff[1]=1;
+      break;
+    }
+    default:{
+      par_print_err(Myrank,"unknown type of aggregation is required",__FILE__,__LINE__,__func__);
+    }
+    }
+    for (i=1;i<Nproc;i++){
+      if (elem_domain == aggregate_domain)
+        buff[0]=neldom[i];
+      MPI_Send (buff,2,MPI_LONG,i,Myrank,MPI_COMM_WORLD);
+      //MPI_Send_mine (buff,2,MPI_LONG,i,Myrank,MPI_COMM_WORLD);
+    }
+    if (elem_domain == aggregate_domain)
+      buff[0]=neldom[0];
+  }
+  else{
+    MPI_Recv (buff,2,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
+    //MPI_Recv_mine (buff,2,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
+  }
+  MPI_Barrier (MPI_COMM_WORLD);
+  //  the number of connections between the macroproblem and microproblem
+  //  the number of entities (elements or subdomains) of macroproblem connected with microproblem
+  nelidom=buff[0];
+  //  the maximum number of connections between the macroproblem and microproblem
+  //  the maximum number of entities (elements or subdomains) of macroproblem connected with microproblem
+  buffsize1=buff[1];
+  //  dimension of the problem, one, two or three dimensional case
+  //  the dimension is defined with respect to the first finite element
+  dim = Mt->give_dimension (0);
+  //  the number of stress/strain components
+  ncomp = Mt->give_ncomp (0,0);
+  buffsize2=buffsize1;
+  //  buffer for displacments and strain components
+  buffsize1*=(dim+ncomp+1);
+  //  buffer for the stiffness matrix
+  buffsize2*=(ncomp*ncomp);
+  fprintf (stdout,"\n procesor %d  nelidom   %ld",Myrank,nelidom);
+  fprintf (stdout,"\n procesor %d  buffsize1 %ld",Myrank,buffsize1);
+  fprintf (stdout,"\n procesor %d  buffsize2 %ld",Myrank,buffsize2);
+  delete [] buff;
+  buff1 = new double [buffsize1];
+  nullv (buff1,buffsize1);
+  buff2 = new double [buffsize2];
+  nullv (buff2,buffsize2);
+  /**************************************************************************/
+  if (Myrank==0){
+    //  the number of unknonws
+    //  it is the number of unknowns in macroproblem in the case of the master processor
+    //  otherwise, it is the number of unknowns in microproblems
+    //  number of mechanical degrees of freedom
+    n=Ndofm;
+    //  vector of prescribed force loads, it does not contain forces caused by temperature, etc.
+    fl   = new double [n];
+    nullv (fl,n);
+    //  vector of resulting internal forces
+    fi   = new double [n];
+    nullv (fi,n);
+    rhs = Lsrs->give_rhs (0);
+  }
+  // fictious time
+  end_time = 0.0;
+  if (Myrank==0){
+    for (ii=1;ii<Nproc;ii++){
+      MPI_Send (&end_time,1,MPI_DOUBLE,ii,Myrank,MPI_COMM_WORLD);
+      //MPI_Send_mine (&end_time,1,MPI_DOUBLE,ii,Myrank,MPI_COMM_WORLD);
+    }
+  }
+  else{
+    MPI_Recv (&end_time,1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
+    //MPI_Recv_mine (&end_time,1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
+  }
+  MPI_Barrier (MPI_COMM_WORLD);
+  fprintf (stdout,"\n procesor %d",Myrank);
+  if (Myrank==0){
+    // **********************************
+    //  solution of microscale problems
+    // **********************************
+    for (ii=1;ii<Nproc;ii++){
+      //  function assembles data from the macroproblem for the ii-th microproblem to the array buff1
+      // tady pole buff1 musi mit rozmer podle nelidom??!! nebo podle poctu elementu??!!
+      higher_to_lower_levelm (ii,buffsize1,buff1);
+      MPI_Send (buff1,buffsize1,MPI_DOUBLE,ii,Myrank,MPI_COMM_WORLD);
+      //MPI_Send_mine (buff1,buffsize1,MPI_DOUBLE,ii,Myrank,MPI_COMM_WORLD);
+    }
+  }
+  else{
+    MPI_Recv (buff1,buffsize1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
+    //MPI_Recv_mine(buff1,buffsize1,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
+    for (ii=0;ii<nelidom;ii++){
+      fprintf (Out,"\n\n *********************************************");
+      fprintf (Out,"\n Myrank  = %ld",Myrank);
+      fprintf (Out,"\n nelidom = %ld",ii);
+      fprintf (Out,"\n\n pred vypoctem paral_mechanical_homogenization:");
+      for (long ijk=0;ijk<buffsize1;ijk++){
+        //for (long ijk=0;ijk<1;ijk++){
+          fprintf (Out,"\n buff1 %ld   %le",ijk,buff1[ijk]);
+      }
+      fprintf (Out,"\n");
+      for (long ijk=0;ijk<buffsize2;ijk++){
+        //for (long ijk=0;ijk<1;ijk++){
+        fprintf (Out,"\n buff2 %ld   %le",ijk,buff2[ijk]);
+      }
+      fflush(Out);
+      paral_mechanical_homogenization (buff1+ii*(dim+ncomp+1),buff2+ii*(ncomp*ncomp));
+      fprintf (Out,"\n\n po vypoctu paral_mechanical_homogenization:");
+      for (long ijk=0;ijk<buffsize1;ijk++){
+        //for (long ijk=0;ijk<1;ijk++){
+        fprintf (Out,"\n buff1 %ld   %le",ijk,buff1[ijk]);
+      }
+      fprintf (Out,"\n");
+      for (long ijk=0;ijk<buffsize2;ijk++){
+        //for (long ijk=0;ijk<1;ijk++){
+        fprintf (Out,"\n buff2 %ld   %le",ijk,buff2[ijk]);
+      }
+      fflush(Out);
+    }
+  }
+  MPI_Barrier (MPI_COMM_WORLD);
+  /**************************************************************************/
+  if (Myrank==0){
+    for (ii=1;ii<Nproc;ii++){
+      MPI_Recv (buff2,buffsize2,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
+      //MPI_Recv_mine(buff2,buffsize2,MPI_DOUBLE,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
+      k=stat.MPI_TAG;
+      jj=0;
+      for (j=0;j<Mt->ne;j++){
+        //fprintf (Out,"\n Myrank  = %ld",Myrank);
+        //fprintf (Out,"\n element no. = %ld",j);
+        //fprintf (Out,"\n Gtm->stop->eldom[j] = %ld",Gtm->stop->eldom[j]);
+        //fprintf (Out,"\n ii = %ld",ii);
+        //fprintf (Out,"\n k  =stat.MPI_TAG = %ld",k);
+        if (Gtm->stop->eldom[j]==k){
+          //  only elements connected with the k-th microscale are assumed
+          switch (Mp->mami){
+          case elem_domain:{
+            //  each element is connected with its own microproblem
+            Mm->hommatm[j].assemble_matrices (buff2,ncomp,dim);
+            break;
+          }
+          case aggregate_domain:{
+            //  generation of auxiliary data needed in parallel computation
+            Mm->hommatm[j].assemble_matrices (buff2+jj*(ncomp*ncomp),ncomp,dim);
+            jj++;
+            break;
+          }
+          case material_aggregate_domain:{
+            //  aggregate sends averaged data to microproblem
+            //  index jj is not incremented in this case because all elements in aggregate??!!
+            //  obtain the same conductivity and capacity matrices
+            //
+            //  index of subdomain-aggregate = index of material (processor)
+            //fprintf (Out,"\n Myrank  = %ld",Myrank);
+            //fprintf (Out,"\n from Myrank  = %ld",ii);
+            //fprintf (Out,"\n element  = %ld",j);
+            //fprintf (Out,"\n\n pred ulozenim matice na materialu c. %ld:",j);
+            //for (long ijk=0;ijk<buffsize2;ijk++){
+            //fprintf (Out,"\n buff2 %ld   %le",ijk,buff2[ijk]);
+            //}
+            Mm->hommatm[j].assemble_matrices (buff2,ncomp,dim);
+            //fflush(Out);
+            break;
+          }
+          case eff_aggregate_domain:{
+            //  aggregate sends averaged data to microproblem
+            //  index jj is not incremented in this case because all elements in aggregate
+            //  obtain the same conductivity and capacity matrices
+            //
+            //  index of subdomain-aggregate = index of element material
+            kk = Mt->elements[i].idm[0];
+            Mm->hommatm[kk].assemble_matrices (buff2+kk*(ncomp*ncomp),ncomp,dim);
+            break;
+          }
+          default:{
+            par_print_err(Myrank,"unknown type of aggregation is required",__FILE__,__LINE__,__func__);
+          }
+          }
+        }
+      }//  end loop over the number of elements
+    }//  end loop over the number of processors
+  }//  end of the master processor part
+  else{
+    MPI_Send (buff2,buffsize2,MPI_DOUBLE,0,Myrank,MPI_COMM_WORLD);
+    //MPI_Send_mine (buff2,buffsize2,MPI_DOUBLE,0,Myrank,MPI_COMM_WORLD);
+  }
+  MPI_Barrier (MPI_COMM_WORLD);
+  fflush (stdout);
+  if (Myrank==0){
+    //  stiffness matrix assembling
+    stiffness_matrix (0);
+    //  loop over load cases
+    for (i=0;i<Lsrs->nlc;i++){
+      //  right hand side of system
+      mefel_right_hand_side (i,rhs);
+    }
+    Smat->printmat (Out);
+    //  solution of equation system
+    for (i=0;i<Lsrs->nlc;i++){
+      Mp->ssle->solve_system (Gtm,Smat,Lsrs->give_lhs(i),Lsrs->give_rhs(i),Out);
+    }
+    print_init(-1, "wt");
+    // clean temperature strains because they are stored in one array used for all load cases
+    Mm->nulltempstrains();
+    rhs = Lsrs->give_rhs(0);
+    for (i=0;i<Lsrs->nlc;i++){
+      //  indicator of strain computation
+      //  it means, strains at integration points have not been computed
+      Mp->strainstate=0;
+      //  indicator of stress computation
+      //  it means, stresses at integration points have not been computed
+      Mp->stressstate=0;
+      //  indicator of computation of other array
+      //  it means, stresses at integration points have not been computed
+      Mp->otherstate=0;
+      // left hand side vector for debugging
+      lhs =  Lsrs->give_lhs(i);
+      // rhs was rewritten by the solver call -> rebuild load vector
+      aux = Mp->homog;
+      Mp->homog = 0;
+      mefel_right_hand_side(i,rhs,fl);
+      Mp->homog = aux;
+      //  computes required quantities
+      compute_req_val (i);
+      // in artificial time 0.0 the nodal forces contains only the load vector fl
+      print_step(i, 0, 0.0, fl);
+      // if the nodal forces are required for graphical output
+      // reaction + load vector are printed in the artificial time 1.0
+      if (Outdm->nog.selnforce.st != sel_no)
+        {
+          //  indicator of strain computation
+          //  it means, strains at integration points have not been computed
+          Mp->strainstate=0;
+          //  indicator of stress computation
+          //  it means, stresses at integration points have not been computed
+          Mp->stressstate=0;
+          //  indicator of computation of other array
+          //  it means, stresses at integration points have not been computed
+          Mp->otherstate=0;
+          // calculate nodal forces and reactions
+          internal_forces(i, fi);
+          //  computes required quantities
+          compute_req_val (i);
+          //  prints required quantities
+          // in artificial time 1.0 the nodal forces contains the load vector + reactions
+          Outdm->print_graphics(Outdm->outgr, i, 1.0, 1, fi);
+        }
+      // clean temperature strains because they are stored in one array used for all load cases
+      Mm->nulltempstrains();
+    }
+  }
+  fprintf (stdout,"\n procesor %d",Myrank);
+  fflush (stdout);
+  if (Myrank==0){
+    print_close();
+    delete [] fl;
+    delete [] fi;
+  }
+  delete [] buff1;
+  delete [] buff2;
+}
+/**
+   function selects appropriate components of higher level of a multiscale
+   problem which are sent to lower level
+   @param llid - lower level id
+   @param ncbuff - the number of components in the array buff
+   @param buff - array for selected components
+   TKr 14/07/2022 according to JK
+*/
+void higher_to_lower_levelm (long llid,long ncbuff,double *buff)
+{
+  long i,j,nc,*counter,ncomp,dim;
+  double area,ar,*aux;
+  strastrestate ssst;
+  counter = new long [1];
+  counter[0]=0;
+  switch (Mp->mami){
+  case elem_domain:{
+    //  each element is connected with a microproblem
+    //  loop over the number of elements
+    for (i=0;i<Mt->ne;i++){
+      if (Gtm->stop->eldom[i]==llid){
+        //  only elements connected with the microid-th microscale are assumed
+        //  the following function is located in the file MEFEL/SRC/elemswitch
+        higher_to_lower_level_elemm (i,counter,buff);
+      }
+    }//  end of loop over the number of elements
+    break;
+  }
+  case aggregate_domain:{
+    //  each element is connected with a microproblem
+    //  elements in one aggregate are connected with the same microproblem
+    //  loop over the number of elements
+    for (i=0;i<Mt->ne;i++){
+      if (Gtt->stop->eldom[i]==llid){
+        //  the following function is located in the file MEFEL/SRC/elemswitch
+        higher_to_lower_level_elemm (i,counter,buff);
+      }
+    }//  end of loop over the number of elements
+    break;
+  }
+  case material_aggregate_domain:{
+    //  each aggregate (elements of one material) is connected with a microproblem
+    //  averaged data is sent to the microproblem
+    //  dimension of the problem, one, two or three dimensional case
+    //  the dimension is defined with respect to the first finite element
+    dim = Mt->give_dimension (0);
+    //  the number of transported materials
+    ncomp = Mt->give_ncomp (0,0);
+    aux = new double [ncbuff];
+    for (i=0;i<ncbuff;i++){
+      buff[i]=0.0;//erase actual position of array buff
+    }
+    // total area of the aggregate
+    ar = 0.0;
+    //  loop over the number of elements
+    for (i=0;i<Mt->ne;i++){
+      //  generally, the averaging should take into account the size
+      //  of particular finite elements, not only their number
+      //  if needed, the function higher_to_lower_level_elem
+      //  will return also area or volume of finite elements
+      if (Gtm->stop->eldom[i]==llid){
+        //  the following function is located in the file TRFEL/SRC/elemswitcht
+        counter[0]=0;
+        for (long ijk=0;ijk<ncbuff;ijk++){
+          aux[ijk]=0.0;//erase array
+        }
+        higher_to_lower_level_elemm (i,counter,aux);//array aux is filled with displacements and their strains
+        //averaging according to element area/volume:
+        ssst=Mt->give_ssst (i,0);
+        if(ssst == spacestress)//for 3D elements
+          area = Mt->give_volume(i);
+        else
+          area = Mt->give_area(i);
+        area = fabs(area);
+        //fprintf (Out,"\n elem %ld area = %lf",i,area);
+        for (j=0;j<ncbuff;j++){
+          buff[j]+=aux[j]*area;
+        }
+        ar += area;
+      }
+    }//  end of loop over the number of elements
+    //fprintf (Out,"\n\n pred prumerovanim buff1:");
+    //fprintf (Out,"\n Proc No.= llid = %ld po: \n",llid);
+    //for (long ijk=0;ijk<ncbuff;ijk++){
+    //  fprintf (Out,"\n buff1 %ld   %le",ijk,buff[ijk]);
+    //}
+    //  averaging of the values
+    fprintf (stdout,"\n\n Averaging of values of material aggregate on processor No. %ld",llid);
+    for (i=0;i<ncbuff;i++){
+      buff[i]=buff[i]/ar;
+    }
+    //fprintf (Out,"\n\n zprumerovane buff1:");
+    //fprintf (Out,"\n Proc No.= llid = %ld po: \n",llid);
+    //for (long ijk=0;ijk<ncbuff;ijk++){
+    //fprintf (Out,"\n buff1 %ld   %le",ijk,buff[ijk]);
+    //}
+    delete [] aux;
+    break;
+  }
+  case eff_aggregate_domain:{
+    //  each aggregate is connected with a microproblem
+    //  averaged data is sent to the microproblem
+    aux = new double [ncbuff];
+    for (i=0;i<ncbuff;i++){
+      aux[i]=0.0;
+      buff[i]=0.0;
+    }
+    //  the number of contributions = the number of elements in the aggregate
+    nc=0;
+    //  loop over the number of elements
+    for (i=0;i<Mt->ne;i++){
+      //  generally, the averaging should take into account the size
+      //  of particular finite elements, not only their number
+      //  if needed, the function higher_to_lower_level_elem
+      //  will return also area or volume of finite elements
+      if (Gtm->stop->eldom[i]==llid){
+        //  the following function is located in the file MEFEL/SRC/elemswitch
+        higher_to_lower_level_elemm (i,counter,buff);
+        counter[0]=0;
+        //nc++;
+        //averaging according to element area/volume:
+        ssst=Mt->give_ssst (i,0);
+        if(ssst == spacestress)//for 3D elements
+          area = Mt->give_volume(i);
+        else
+          area = Mt->give_area(i);
+        area = fabs(area);
+        for (j=0;j<ncbuff;j++){
+          aux[j]+=buff[j]/area;
+        }
+      }
+    }//  end of loop over the number of elements
+    //  averaging of the values
+    //for (i=0;i<nc;i++){
+    //buff[i]=aux[i]/nc;
+    //}
+    delete [] aux;
+    break;
+  }
+  default:{
+    par_print_err (Myrank,"unknown type of aggregation is required",__FILE__,__LINE__,__func__);
+  }
+  }
+  delete [] counter;
+}

 /trunk/SIFEL/HPARAL/SRC/hplssolver.h
 ,6 → 3,8
 //#include <stdio.h>
+void parallel_homogenization_linear_statics_tiles ();
 void parallel_homogenization_linear_statics ();
+void higher_to_lower_levelm (long llid,long ncbuff,double *buff);
 #endif

 /trunk/SIFEL/HPARAL/SRC/hpmefel.cpp
 ,8 → 12,6
 int main (int argc,char *argv[])
 {
   time_t bt,et;
-  stochdriver *stochd;
-  stochd = new stochdriver;
   bt = time (NULL);
 ,7 → 22,8
   MPI_Get_processor_name(proc_namet,&nameLengtht);
   //  initialization of the code
-  hpmefel_init (argc,(const char **)(argv),stochd);
+  hpmefel_init (argc,(const char **)(argv));
+  //hpmefel_init_tiles (argc,(const char **)(argv));
   //  solution of the problem
   solve_mefel_problem_parallel ();
 ,17 → 30,17
   et = time (NULL);
-  /*
-  fprintf (Out,"\n\n\n Udaje o dobach vypoctu \n");
-  fprintf (Out,"\n\n celkova doba vypoctu             %ld",et-bt);
+  //fprintf (Out,"\n\n\n Udaje o dobach vypoctu \n");
+  //fprintf (Out,"\n\n celkova doba vypoctu             %ld",et-bt);
   if (Myrank==0){
     fprintf (stdout,"\n\n\n Udaje o dobach vypoctu \n");
     fprintf (stdout,"\n\n celkova doba vypoctu             %ld",et-bt);
   }
-  fprintf (Out,"\n");
-  fclose (Out);
-  */
+  //fprintf (Out,"\n");
+  //fclose (Out);
   fprintf (stdout,"\n");  fprintf (stderr,"\n");
   MPI_Finalize ();

 /trunk/SIFEL/HPARAL/SRC/hpmefelinit.cpp
 ,7 → 1,6
 #include "hpmefelinit.h"
 #include "hpglobal.h"
 #include "xfile.h"
-#include "stochdriver.h"
 #include "stacktrace.h"
 /**
 ,11 → 7,69
    function reads all data from the input file
    @param argv - pointer to input file names
-   @param stochd - pointer to stochdriver
+   TKr, 26/07/2022
+*/
+void hpmefel_init (int argc,const char *argv[])
+{
+  long i,mp;
+  char name[1100];
+  MPI_Status stat;
+  const char *aux;
+  //  names of input files
+  strcpy (name,argv[1]);
+  sprintf (&name[strlen (argv[1])],"%d.in",Myrank+1);
+  aux = argv[1];
+  argv[1] = name;
+  //  reading of the classical MEFEL input file
+  mefel_init (argc,(const char **)(argv));
+  argv[1] = aux;
+  if (Myrank==0){
+    switch (Mp->tprob){
+    case linear_statics:{
+      mp=1;
+      break;
+    }
+    case mech_timedependent_prob:{
+      mp=15;
+      break;
+    }
+    default:{
+      print_err("unknown type of mechanical problem is required",__FILE__,__LINE__,__func__);
+    }
+    }
+    for (i=1;i<Nproc;i++){
+      MPI_Send (&mp,1,MPI_LONG,i,Myrank,MPI_COMM_WORLD);
+    }
+  }
+  else{
+    MPI_Recv (&mp,1,MPI_LONG,MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);
+  }
+  MPI_Barrier (MPI_COMM_WORLD);
+  if (mp==1){
+    Mtprobm=linear_statics;
+  }
+  if (mp==15){
+    Mtprobm=mech_timedependent_prob;
+  }
+}
+/**
+   function reads all data from the input file
+   @param argv - pointer to input file names
    JK, 4.6.2013
+   TKr, actualized 12/07/2022
 */
-void hpmefel_init (int argc,const char *argv[], stochdriver *stochd)
+void hpmefel_init_tiles (int argc,const char *argv[])
 {
   long i,j,lcid;
   char name[1100],namebackup[1100];
 ,7 → 179,7
     argv[1] = name;
     //  reading of the classical MEFEL input file
-    mefel_init (argc,(const char **)(argv),stochd);
+    mefel_init (argc,(const char **)(argv));
     Nbdof[i]=Ndofm-Nidof[i];
 ,6 → 224,10
       Mtprobm=mat_nonlinear_statics;
       break;
     }
+    case mech_timedependent_prob:{
+      Mtprobm=mech_timedependent_prob;
+      break;
+    }
     default:{
     }

 /trunk/SIFEL/HPARAL/SRC/hpmefelinit.h
 ,8 → 1,8
 #ifndef HPMEFELINIT_H
 #define HPMEFELINIT_H
-class stochdriver;
-void hpmefel_init (int argc,const char *argv[], stochdriver *stochd);
+void hpmefel_init (int argc,const char *argv[]);
+void hpmefel_init_tiles (int argc,const char *argv[]);
 #endif

 /trunk/SIFEL/HPARAL/SRC/hpnpsolvert.cpp
 ,7 → 238,6
     print_stept(0,i,Tp->time,NULL);
     print_flusht();
   }
   /**************************************************************************/
   if (Myrank==0){
 ,7 → 253,7
   MPI_Barrier (MPI_COMM_WORLD);
-  fprintf (stdout,"\n procesor %d end_time %lf",Myrank,end_time);
+  fprintf (stdout,"\n procesor %d end_time %lf\n",Myrank,end_time);
   do{
 ,8 → 492,6
     }
     fprintf (stdout,"\n procesor %d time = %lf end_time = %lf",Myrank,time,end_time);
     fflush (stdout);
-    MPI_Barrier (MPI_COMM_WORLD);
   }while(time<end_time);

 /trunk/SIFEL/HPARAL/SRC/hpsolverm.cpp
 ,5 → 1,6
 #include "hpsolverm.h"
 #include "hplssolver.h"
+#include "hpmtsolver.h"
 #include "hpglobal.h"
 #include "global.h"
 #include "iotools.h"
 ,6 → 9,7
 {
   switch (Mtprobm){
   case linear_statics:{
+    //parallel_homogenization_linear_statics_tiles ();
     parallel_homogenization_linear_statics ();
     break;
   }
 ,6 → 16,10
   case mat_nonlinear_statics:{
     break;
   }
+  case mech_timedependent_prob:{
+    parallel_homogenization_solve_time_dep_prob ();
+    break;
+  }
   default:{
     par_print_err(Myrank,proc_namet,"unknown problem type is required",__FILE__,__LINE__,__func__);
   }

/trunk/SIFEL/HPARAL/SRC/hpsolverm.h
2,5 → 2,6
#define HPSOLVERM_H

void solve_mefel_problem_parallel ();
long compare_on_master (double val, double lim);

#endif

 /trunk/SIFEL/HPARAL/SRC/hptrfel.cpp
 ,8 → 12,6
 int main (int argc,char *argv[])
 {
   time_t bt,et;
-  stochdrivert *stochd;
-  stochd = new stochdrivert;
   XFILE *in;
   bt = time (NULL);
 ,7 → 23,7
   MPI_Get_processor_name(proc_namet,&nameLengtht);
   //  initialization of the code
-  hptrfel_init (argc,(const char **)(argv),stochd);
+  hptrfel_init (argc,(const char **)(argv));
   //hptrfel_init_tiles (argc,(const char **)(argv),stochd);
   //hptrfel_init_tiles_parsolver (argc,(const char **)(argv),stochd, in);

 /trunk/SIFEL/HPARAL/SRC/hptrfelinit.cpp
 ,7 → 1,6
 #include "hptrfelinit.h"
 #include "hpglobal.h"
 #include "xfile.h"
-#include "stochdrivert.h"
 #include "stacktrace.h"
 /**
 ,11 → 7,11
    function reads all data from the input file
    @param argv - pointer to input file names
-   @param stochd - pointer to stochdriver
    JK, 22.9.2010
+   TKr, actualized 12/07/2022
 */
-void hptrfel_init (int argc,const char *argv[], stochdrivert *stochd)
+void hptrfel_init (int argc,const char *argv[])
 {
   long i,mpt;
   char name[1100];
 ,7 → 25,7
   argv[1] = name;
   //  reading of the classical TRFEL input file
-  trfel_init (argc,(const char **)(argv),stochd);
+  trfel_init (argc,(const char **)(argv));
   argv[1] = aux;
 ,7 → 44,7
       break;
     }
     default:{
+      print_err("unknown type of transport problem is required ",__FILE__, __LINE__, __func__);
     }
     }
 ,11 → 73,11
    function reads all data from the input file
    @param argv - pointer to input file names
-   @param stochd - pointer to stochdriver
    JK, 4.6.2013
+   TKr, actualized 12/07/2022
 */
-void hptrfel_init_tiles (int argc,const char *argv[], stochdrivert *stochd)
+void hptrfel_init_tiles (int argc,const char *argv[])
 {
   long i,j,lcid;
   char name[1100],namebackup[1100];
 ,7 → 216,7
     argv[1] = name;
     //  reading of the classical MEFEL input file
-    trfel_init (argc,(const char **)(argv),stochd);
+    trfel_init (argc,(const char **)(argv));
     Nbdof[i]=Ndoft-Nidof[i];
 ,11 → 269,11
    function reads all data from the input file
    @param argv - pointer to input file names
-   @param stochd - pointer to stochdriver
    JK, 9. 9. 2014
+   TKr, actualized 12/07/2022
 */
-void hptrfel_init_tiles_parsolver (int argc,const char *argv[], stochdrivert *stochd, XFILE* &in)
+void hptrfel_init_tiles_parsolver (int argc,const char *argv[], XFILE* &in)
 {
   long i,lcid;
   char name[1100];
 ,7 → 302,7
   argv[1] = name;
   //  reading of the classical TRFEL input file
-  trfel_init (argc,(const char **)(argv),stochd);
+  trfel_init (argc,(const char **)(argv));
   //  restoration of the original name in argv[1]
   argv[1]=bargv;

 /trunk/SIFEL/HPARAL/SRC/hptrfelinit.h
 ,13 → 1,12
 #ifndef HPTRFELINIT_H
 #define HPTRFELINIT_H
-class stochdrivert;
 struct XFILE;
-void hptrfel_init (int argc,const char *argv[], stochdrivert *stochd);
+void hptrfel_init (int argc,const char *argv[]);
-void hptrfel_init_tiles (int argc,const char *argv[], stochdrivert *stochd);
+void hptrfel_init_tiles (int argc,const char *argv[]);
-void hptrfel_init_tiles_parsolver (int argc,const char *argv[], stochdrivert *stochd, XFILE* &in);
+void hptrfel_init_tiles_parsolver (int argc,const char *argv[], XFILE* &in);
 #endif

 /trunk/SIFEL/MEFEL/SRC/Makefile
 ,7 → 41,7
 SRCS +=  elastisomat.cpp elastisopdmat.cpp elastortomat.cpp elasttime.cpp element.cpp elemparticle.cpp elemswitch.cpp endnodem.cpp epsolver.cpp fdsolver.cpp flsubdom.cpp fixortodam.cpp
 SRCS +=  generalmod.cpp
 SRCS +=  geoelast.cpp glasgmech.cpp glasgowdam.cpp global.cpp globmat.cpp
-SRCS +=  graphmat.cpp hardsoft.cpp hissplas.cpp homog.cpp homogmech.cpp hvisolver.cpp inicd.cpp intpoints.cpp isoviscous.cpp hypoplunsatexptherm2.cpp
+SRCS +=  graphmat.cpp hardsoft.cpp hissplas.cpp homog.cpp homogmatm.cpp homogmech.cpp hvisolver.cpp inicd.cpp intpoints.cpp isoviscous.cpp hypoplunsatexptherm2.cpp
 SRCS +=  j2flow.cpp layplate.cpp lbsolver.cpp lemaitre.cpp lenjonesmat.cpp lfssolver.cpp lhsrhs.cpp
 SRCS +=  linhex.cpp linhexrot.cpp lintet.cpp lintetrot.cpp linwedge.cpp llssolver.cpp loadcase.cpp

 /trunk/SIFEL/MEFEL/SRC/alias.h
 ,7 → 148,7
 // aliases for material models
-enum mattype {elisomat=1,elgmat3d=2,elortomat=3,elgmat2d=4,
+enum mattype {elisomat=1,elgmat3d=2,elortomat=3,elgmat2d=4, homomatm=9,
               simplas1d=10,jflow=11,mohrcoulparab=22,mohrcoul=24,
               boermaterial=25,druckerprager=26,doubledrprager=27,druckerprager2=28,modcamclaymat=30, modcamclaycoupmat=31,bbmcoupmat=33,
               shefpl=40,chenplast=42, hissplasmat=45,
 ,7 → 167,7
               damage_plasticity=1000, viscoplasticity=1010, viscoelasticity=1020,
               creep_damage=1030, time_switchmat=1040, effective_stress=1050, shrinkagemat=1060, elasttimemat=1070,
               lenjonespot=2000, cusatismat=2010};
-const enumstr mattypestr[] = {{"elisomat",1}, {"elgmat3d",2}, {"elortomat",3}, {"elgmat2d",4},{"simplas1d",10},
+const enumstr mattypestr[] = {{"elisomat",1}, {"elgmat3d",2}, {"elortomat",3}, {"elgmat2d",4},{"homomatm",9},{"simplas1d",10},
                               {"jflow",11},  {"mohrcoulparab",22},
                               {"mohrcoul",24}, {"boermaterial",25}, {"druckerprager",26},  {"doubledrprager",27}, {"druckerprager2",28},
                               {"modcamclaymat",30}, {"modcamclaycoupmat",31}, {"bbmcoupmat",33}, {"shefpl",40}, {"chenplast",42}, {"hissplasmat",45},

 /trunk/SIFEL/MEFEL/SRC/elastisomat.cpp
 ,6 → 394,7
 */
 void elastisomat::matstiff_spacestr (matrix &d)
 {
+  long i,j;
   double g,s;
   fillm(0.0,d);
 ,6 → 407,16
   d[2][0]=d[0][1];   d[2][1]=d[0][1];  d[2][2]=d[0][0];
   d[3][3]=g;         d[4][4]=g;        d[5][5]=g;
+  ////////////////////////////////////////////////////////////
+  // this printing below is only for debug:
+  fprintf(Out,"\n\n vysledna matice D^M\n");
+  for (i=0;i<6;i++){
+    for (j=0;j<6;j++)
+      fprintf(Out,"%e  \n",d[i][j]);
+    fprintf(Out,"\n");
+  }
 }

 /trunk/SIFEL/MEFEL/SRC/elemswitch.cpp
 ,3 → 5723,58
     }
   }
 }
+/**
+   function selects components from the macro level to meso level for homogenization
+   @param eid - element id
+   @param counter - actual position in the array buff
+   @param buff - array containing components
+   TKr 22/07/2022 accroding to JK
+*/
+void higher_to_lower_level_elemm (long eid,long *counter,double *buff)
+{
+  elemtype te;
+  //  type of element
+  te=Mt->give_elem_type (eid);
+  switch (te){
+    case planeelementqt:
+      //Peqt->higher_to_lower_level (eid,counter,buff);
+      break;
+    case planeelementrotlt:
+      //Perlt->higher_to_lower_level (eid,counter,buff);
+      break;
+    case planeelementlq:
+      //Pelq->higher_to_lower_level (eid,counter,buff);
+      break;
+    case planeelementqq:
+      //Peqq->higher_to_lower_level (eid,counter,buff);
+      break;
+    case axisymmlq:
+      //Asymlq->higher_to_lower_level (eid,counter,buff);
+      break;
+    case axisymmqq:
+      //Asymqq->higher_to_lower_level (eid,counter,buff);
+      break;
+    case lineartet:
+      //Ltet->higher_to_lower_level (eid,counter,buff);
+      break;
+    case quadrtet:
+      //Qtet->higher_to_lower_level (eid,counter,buff);
+      break;
+    case linearhex:
+      //Lhex->higher_to_lower_level (eid,counter,buff);
+      break;
+    case quadrhex:
+      //Qhex->higher_to_lower_level (eid,counter,buff);
+      break;
+    default:
+      print_err("unknown element type %d is required", __FILE__, __LINE__, __func__, Mt->give_elem_type(eid));
+  }
+}

 /trunk/SIFEL/MEFEL/SRC/elemswitch.h
 ,4 → 111,8
 /// function computes stresses and eventual state variables at auxiliary integration points
 void compute_aipstresses(long lcid);
+/// function selects components from the macro level to meso level for homogenization problems
+void higher_to_lower_level_elemm (long eid,long *counter,double *buff);
 #endif

 /trunk/SIFEL/MEFEL/SRC/generalmod.cpp
 ,8 → 51,7
   qstatev[6] = SrM
   qstatev[7] = a_scan (required on input)
   qstatev[8] = re
-  qstatev[9] = empty
   Order of parameters in the array rkf_parms
   rkf_parms[0] = err_sig
   rkf_parms[1] = h_min
 ,7 → 76,7
 */
 long Hypoplasti_unsat_expansive_thermal::soil_model(double strain_gen[], double stress_gen[],
-  double qstatev[], double dstrain_gen[], double dtime, double *DDtan_gen, double parms[], double rkf_statev[], double rkf_parms[], int flag, int kinc/*=0*/) {
+  double qstatev[], double dstrain_gen[], double dtime, double *DDtan_gen, double parms[], double rkf_statev[], double rkf_parms[], int flag, int kinc/*=0*/, int ipp/*=0*/) {
   long errorrkf = 0;
   long errormath = 0;
 ,8 → 575,7
   dqstatev[6]=SrM_new-SrM;
   dqstatev[7]=a_scan_new-qstatev[7];
   dqstatev[8]=0; //re: to be filled later
-  dqstatev[9]=0;
   //calculate new state variables, needed just for re
   double qstatev_new[c_nstatev];
   garray_add(qstatev, dqstatev, qstatev_new, c_nstatev);

 /trunk/SIFEL/MEFEL/SRC/generalmod.h
 ,7 → 54,7
   //Updated soil_model with Runge-Kutta schemes
   virtual long soil_model(double strain_gen[], double stress_gen[], double qstatev[],
-        double dstrain_gen[], double dtime, double *DDtan_gen, double parms[], double rkf_statev[], double rkf_parms[], int flag, int kinc=0) {return(0);};
+        double dstrain_gen[], double dtime, double *DDtan_gen, double parms[], double rkf_statev[], double rkf_parms[], int flag, int kinc=0, int ipp=0) {return(0);};
   virtual bool calc_fsigq(double signet[9], double suction, double Temper, double *qstatev,
         double deps[9], double dsuction, double dTemper, double dtime,
 ,7 → 121,7
   static const int c_ngstrain = 8;
   static const int c_ngstress = 7;
-  static const int c_nstatev  = 10;
+  static const int c_nstatev  = 9;
   static const int c_nparms   = 23;
   static const bool debug     = false;
 ,7 → 135,7
   //Updated soil_model with Runge-Kutta schemes
   long soil_model(double strain_gen[], double stress_gen[], double qstatev[],
-        double dstrain_gen[], double dtime, double *DDtan_gen, double parms[], double rkf_statev[], double rkf_parms[], int flag, int kinc=0);
+        double dstrain_gen[], double dtime, double *DDtan_gen, double parms[], double rkf_statev[], double rkf_parms[], int flag, int kinc=0, int ipp=0);
   /* functions which distinguish original and updated versions */
   virtual double give_rlambda(double suction, double a_scan, double dsuction, double sewM, double SrM, double& rlambda_for_ascan, double& rlambda_for_se);

 /trunk/SIFEL/MEFEL/SRC/homogmech.cpp
 ,7 → 25,7
 {
   long i,j,k,ii,jj,hncomp,tnipel,ipp;
-  //  initiation of transport material models
+  //  initiation of mechanical material models
   //Mm->initmaterialmodels();
   //  nodes - integration points interpolation
   //approximation_puc (); //if it is needed
 ,9 → 325,186
                      (in master node on PUC) stored as an array (one dimensional)
    @param matrix_k - homogenized stiffness matrix of PUC stored as an array (one dimensional)
-   TKr, 07/01/2015
+   TKr, 12/07/2022
 */
-void paral_transport_homogenization (double *unkn_r,double *matrix_k)
+void paral_mechanical_homogenization (double *unkn_r,double *matrix_k)
 {
-  //has to be completed??!!
+  long i,j,k,ii,jj,hncomp,tnipel,ipp;
+  //  initiation of mechanical material models
+  //Mm->initmaterialmodels();
+  //  nodes - integration points interpolation
+  //approximation_puc (); //if it is needed
+  //solving of fluctuation parts
+  //only for the first loadcase and elastic material:
+  Mp->stresscomp = 1; //for macrostresses computation
+  Mp->straincomp = 1; //for macrostrains computation
+  switch (Mp->hstrastre){
+  case planestress:{
+    hncomp = 3;
+    break;
+  }
+  case planestrain:{
+    hncomp = 3;
+    break;
+  }
+  case axisymm:{
+    hncomp = 4;
+    break;
+  }
+  case spacestress:{
+    hncomp = 6;
+    break;
+  }
+  default:{
+    print_err("unknown type of strain-stress state for mechanical homogeniaztion is required", __FILE__, __LINE__, __func__);
+  }
+  }//  end of the switch (strastre){
+  //matrix allocation
+  matrix d(ASTCKMAT(hncomp,hncomp));
+  vector mstress(ASTCKVEC(hncomp));
+  //loop over strain-stress components; only for one (first) loadcase
+  for(ii=0;ii<hncomp;ii++){
+    //only macro strain loading:
+    if (Mp->homog == 6)
+      {
+        for(jj=0;jj<hncomp;jj++)//null mstrain
+          Mb->lc[0].mstrain[jj] = 0.0;
+        Mb->lc[0].mstrain[ii] = 1.0;//unit mstrain at ii-th component
+      }
+    else{
+      print_err("only macro strain loading is required",__FILE__,__LINE__,__func__);
+      exit(0);
+    }
+    switch (Mp->tprob){
+    case linear_statics:{
+      solve_linear_statics ();
+      break;
+    }
+    default:{
+      print_err("unknown problem type is required",__FILE__,__LINE__,__func__);
+    }
+    }
+    switch (Mp->homog){
+    case 6:{ // macro-strain approach
+      //stiffness matrix computation from averaged stresses (macrostress)
+      for(i=0;i<hncomp;i++){
+        d[i][ii] = macrostresses(0,i);
+        //fprintf (Out,"macrostresses[%ld][%ld] %16.12le\n",i,ii,d[i][ii]);
+      }
+      //fprintf(Out,"\n Average strain control (%ld):\n",ii+1);
+      for (jj=0; jj<Mm->max_ncompstre; jj++)
+        mstress[jj] = 0.0;
+      for (k=0;k<Mt->ne;k++){
+        if (Gtm->leso[k]==1){
+          if(Mt->give_elem_type (k) == planeelementlq)
+            tnipel = 4;//only the first block for quadrilateral elelemnts
+          else
+            tnipel = Mt->give_tnip(k);
+          ipp = Mt->elements[k].ipp[0][0];
+          for(j=0;j<tnipel;j++){
+            for (jj=0; jj<Mm->max_ncompstre; jj++)
+              mstress[jj]+=Mm->ip[ipp].strain[jj]*Mm->ipv[ipp];
+            //mstress[jj]+=Mm->ipv[ipp];
+            ipp++;
+          }
+        }
+      }
+      for (j=0; j<Mm->max_ncompstre; j++)
+        mstress[j]/=Mt->domvol;
+      //fprintf (Out,"Mt->domvol %16.12le\n",Mt->domvol);
+      //for (jj=0; jj<Mm->max_ncompstre; jj++)
+      //        fprintf(Out,"\n Macrostrain(%ld) = %e\n",jj+1,mstress[jj]);
+      break;
+    }
+    default:{
+      print_err("unknown type of homogenization is required", __FILE__, __LINE__, __func__);
+      abort();
+    }
+    }
+  }//end of loop over strain-stress components
+  //storing of overall properties
+  //Stiffness matrix D
+  k = 0;
+  for (i=0;i<hncomp;i++){
+    for (j=0;j<hncomp;j++){
+      matrix_k[k] = d[i][j];
+      k++;
+    }
+  }
+  ////////////////////////////////////////////////////////////
+  // this printing below is only for debug:
+  /* fprintf(Out,"\n\n vysledna matice D^M\n");
+     for (i=0;i<hncomp;i++){
+     for (j=0;j<hncomp;j++)
+     fprintf(Out,"%e  \n",d[i][j]);
+     fprintf(Out,"\n");
+     }
+  */
+  /*
+    if (Mespr != 0){//printing of matrix into .log file
+    switch (Mp->homog){
+    case 6:{// macro-strain approach
+    fprintf(Out,"\n\n Resulting Stiffness matrix D\n");
+    break;
+    }
+    default:{
+    print_err("unknown type of homogenization is required", __FILE__, __LINE__, __func__);
+    abort();
+    }
+    }
+    for (i=0;i<hncomp;i++){
+    for (j=0;j<hncomp;j++)
+    fprintf(Out,"%e  ",d[i][j]);
+    fprintf(Out,"\n");
+    }
+    //printing for general anisotropic elastic material:
+    switch (Mp->hstrastre){
+    case planestress:
+    case planestrain:{
+    fprintf(Out,"\n\n Print for elastgmat2d:\n");
+    fprintf(Out,"%le %le %le\n",d[0][0],d[0][1],d[0][2]);
+    fprintf(Out,"%le %le\n",d[1][1],d[1][2]);
+    fprintf(Out,"%le\n",d[2][2]);
+    break;
+    }
+    case spacestress:{
+    fprintf(Out,"\n\n Print for elastgmat3d:\n");
+    fprintf(Out,"%le %le %le %le %le %le\n",d[0][0],d[0][1],d[0][2],d[0][3],d[0][4],d[0][5]);
+    fprintf(Out,"%le %le %le %le %le\n",d[1][1],d[1][2],d[1][3],d[1][4],d[1][5]);
+    fprintf(Out,"%le %le %le %le\n",d[2][2],d[2][3],d[2][4],d[2][5]);
+    fprintf(Out,"%le %le %le\n",d[3][3],d[3][4],d[3][5]);
+    fprintf(Out,"%le %le\n",d[4][4],d[4][5]);
+    fprintf(Out,"%le\n",d[5][5]);
+    break;
+    }
+    default:{
+    print_err("unknown type of strain-stress state for mechanical homogeniaztion is required", __FILE__, __LINE__, __func__);
+    }
+    }//  end of the switch (strastre){
+    }
+  */
+  //////////////////////////////////////////////////////////
 }

/trunk/SIFEL/MEFEL/SRC/homogmech.h
4,6 → 4,6
#include <stdio.h>

void mechanical_homogenization();
void paral_mechanical_homogenization(double u,double k,double *c);
void paral_mechanical_homogenization(double unkn_r,double matrix_k);

#endif

 /trunk/SIFEL/MEFEL/SRC/hypoplunsatexptherm2.cpp
 ,7 → 276,7
   eps(6) = statev[1];  // set suction to the corresponding component of generalized strain vector
   eps(7) = statev[3];  // set temperature to the corresponding component of generalized strain vector
   // initialize model
-  herr = huet.soil_model(eps.a, sig.a, statev, deps.a, 1.0, NULL, params, NULL, rkfparams, 0, Mp->jstep);
+  herr = huet.soil_model(eps.a, sig.a, statev, deps.a, 1.0, NULL, params, NULL, rkfparams, 0, Mp->jstep,ipp);
   if (herr)
   {
     print_err("hypoplasticity model cannot be initialized on element %ld (ip=%ld)", __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1);
 ,7 → 290,7
   copyv(statev, Mm->ip[ipp].other+ido+2*ncomp+3, nstatev);
   // set initial values for dSr/ds and dSr/dT according to generalized stiffness matrix
-  herr = huet.soil_model(eps.a, sig.a, statev, NULL, 1.0, dd_gen.a, params, NULL, rkfparams, 3, Mp->jstep);
+  herr = huet.soil_model(eps.a, sig.a, statev, NULL, 1.0, dd_gen.a, params, NULL, rkfparams, 3, Mp->jstep,ipp);
   if (herr)
   {
     print_err("generalized stiffness matrix cannot be initialized on element %ld (ip=%ld)", __FILE__, __LINE__, __func__, Mm->elip[ipp]+1, ipp+1);
 ,7 → 629,10
   rkf_statev[2] = Mm->ip[ipp].eqother[ido+2*ncomp+3+nstatev+7];
   rkf_statev[3] = Mm->ip[ipp].eqother[ido+2*ncomp+3+nstatev];
-  rkf_stat = huet.soil_model(eps.a, sig.a, statev.a, deps.a, dtime, dd_gen.a, params, rkf_statev.a, rkfparams, 1, Mp->jstep);
+  //if((Mp->time >= 5.375400e+05) && ipp == 0)
+  //fprintf (stdout,"\n zde cas!!"); //debug
+  rkf_stat = huet.soil_model(eps.a, sig.a, statev.a, deps.a, dtime, dd_gen.a, params, rkf_statev.a, rkfparams, 1, Mp->jstep, ipp);
   // actualize RKF state variables
   neval = long(rkf_statev[0]);
   nstep = long(rkf_statev[1]);

 /trunk/SIFEL/MEFEL/SRC/lssolver.cpp
 ,7 → 65,13
   fprintf (stdout,"\n\n\n RHS %le\n\n",scprd(rhs,rhs,n));
   */
+  /* fprintf (Out,"\n");
+     fprintf (Out,"\n");
+     for (i=0;i<n;i++){
+     fprintf (Out,"rhs[%ld] %16.12le\n",i,rhs[i]);
+     }
+  */
   //Smat->printmat (Out);
   //  solution of equation system
 ,6 → 79,13
     Mp->ssle->solve_system (Gtm,Smat,Lsrs->give_lhs(i),Lsrs->give_rhs(i),Out);
   }
+  /* fprintf (Out,"\n");
+     fprintf (Out,"\n");
+     for (i=0;i<n;i++){
+     fprintf (Out,"lhs[%ld] %16.12le\n",i,Lsrs->lhs[i]);
+     }
+  */
   /*
   ma=fopen ("reseni.txt","w");
   fprintf (ma,"pocet slozek %ld\n",n);

 /trunk/SIFEL/MEFEL/SRC/mechmat.cpp
 ,6 → 86,7
 #include "inicd.h"
 #include "thervolisomat.h"
 #include "cusatismaterial.h"
+#include "homogmatm.h"
 #include "stacktrace.h"
 #include "ipmap.h"
 #include <math.h>
 ,6 → 163,9
   lplate=NULL;
   cebfipconmat=NULL;
+  //  artificial material obtained from homogenization
+  hommatm=NULL;
   //  material for shrinkage strain computation
   shmat=NULL;
   // damage plasticity interface material
 ,7 → 253,9
   delete [] hypoplustherm;
   delete [] damplifm;
   delete [] cusmat;
+  delete [] hommatm;
   //  material for shrinkage strain computation
   delete [] shmat;
   delete [] eltimemat;
 ,6 → 1127,21
       break;
     }
+    case homomatm:{
+      if (Mespr==1)  fprintf (stdout,"\n  number of material models obtained from homogenization %ld",numt);
+      hommatm = new homogmatm [numt];
+      for (j=0;j<numt;j++){
+        k=numt+1;
+        xfscanf (in,"%ld",&k);
+        if (k>numt || k<1){
+          print_err("wrong number of material obtained from homogenization is required",__FILE__,__LINE__,__func__);
+          abort();
+        }
+        hommatm[k-1].read (in);
+      }
+      break;
+    }
     case simplas1d:{
       if (Mespr==1)  fprintf (stdout,"\n number of simple plastic 1D materials   %ld",numt);
       spl1d = new splas1d [numt];
 ,6 → 2629,9
       eliso[ip[ipp].idm[im]].initval(ipp, im, ido);
       break;
     }
+    case homomatm:
+      hommatm[ip[ipp].idm[im]].initval(ipp, im, ido);
+      break;
     case elortomat:
     case microplaneM4:
     case nonlocalmicroM4:
 ,6 → 3638,8
       ncompo = givencompeqother (ipp,im);
       break;
     }
+    case homomatm:
+      break;
     case creeprs:
     case creepb3:
     case creepdpl:{
 ,6 → 3800,8
       ncompo=ncompstr+13;
       ncompo += givencompeqother (ipp, im+1);
       break;
+    case homomatm:
+      break;
     case hypoplastmat:
       ncompo=2*ncompstr+hypopl[ip[ipp].idm[im]].nstatv+1;
       break;
 ,6 → 5780,7
     case varelisomat:
     case elisopdmat:
     case geoelast:
+    case homomatm:
     case simplas1d:
     case jflow:
     case microplaneM4:
 ,6 → 6090,11
     bbm[i].matstiff (d,ipp,ido);
     break;
   }
+  case homomatm:{
+    i=ip[ipp].idm[im];
+    hommatm[i].matstiff (d,ip[ipp].ssst);
+    break;
+  }
   case hypoplastmat:{
     i=ip[ipp].idm[im];
     hypopl[i].matstiff (d,ipp,ido);
 ,6 → 6673,10
   case bbmcoupmat:{
     bbm[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
     break;
+  }
+  case homomatm:{
+    hommatm[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
+    break;
   }
   case hypoplastmat:{
     hypopl[ip[ipp].idm[im]].nlstresses (ipp,im,ido);
 ,6 → 6926,7
     break;
   }
   case elisomat:
+  case homomatm:
   case simplas1d:
   case jflow:
   case microplaneM4:
 ,6 → 7877,9
   case elisomat :
     eliso[ip[ipp].idm[im]].updateval (ipp,ido);
     break;
+  case homomatm:
+    hommatm[ip[ipp].idm[im]].updateval(ipp,ido);
+    break;
   case elortomat :
     break;
   case winklerpasternak:
 ,9 → 8404,10
   long ncompo;
   switch (ip[ipp].tm[im])
-  {
+    {
     case elisomat :
     case elortomat :
+    case homomatm:
     case simplas1d:
     case jflow:
     case mohrcoul:
 ,6 → 9732,7
   //fprintf (Out,"\n\n int. point %ld \n",ipp);
   //  main iteration loop
   for (i=0;i<ni;i++){
 ,6 → 11712,7
   */
   case elisomat:
   case elortomat:
+  case homomatm:
   case simplas1d:
   case simvisc:
   case jflow:

 /trunk/SIFEL/MEFEL/SRC/mechmat.h
 ,11 → 8,13
 #include "iotools.h"
 #include "selection.h"
 #include "intpoints.h"
+#include "homogmatm.h"
 class elastisomat;
 class elastortomat;
 class elastgmat3d;
 class elastgmat2d;
+class homogmatm;
 class splas1d;
 class j2flow;
 class microM4;
 ,6 → 849,9
   ///  elastic general 2D material
   elastgmat2d *elgm2d;
+  //  artificial material obtained from homogenization
+  homogmatm *hommatm;
   ///  simple plastic 1D material
   splas1d *spl1d;

/trunk/SIFEL/MEFEL/SRC/mtsolver.cpp
810,6 → 810,7

// perform one load step with the help of Newton-Raphson procedure
j = 0;

norf = gnewton_raphson_one_step(lcid, Mp->nlman, fl, r, fb, dr, fi, dtr, istep-1, j, li, fusm);

// just for evaluation of hypoplasticity model ???!!!

 /trunk/SIFEL/MEFEL/SRC/newtonraph.cpp
 ,7 → 706,7
 double gnewton_raphson_one_step(long lcid, nonlinman *nlman, double *fa, double *ra, double *fb, double *dr, double *fi,
                                 double &dtr, long istep, long &j, long li, answertype fusm)
 {
-  long n = Ndofm;    // number of rows of the matrix or vectors
+  long i,n = Ndofm;    // number of rows of the matrix or vectors
   long ini;          // maximum number of inner iterations
   double ierr;       // required error of the residual vector
   double norf;       // norm of residual vector (vector of unbalanced forces)
 ,12 → 724,28
   // type of residual vector norm
   normt = nlman->rnormtnr;
+  //if(Mp->time >= 5.825400e+05)
+  //fprintf (stdout,"\n zde cas!!"); //debug
   //  assembling of tangent stiffness matrix
   assemble_stiffness_matrix(lcid,istep,-1,li,fusm);
+  /* if (Mp->time >= 5.375400e+05 && Mp->time<= 5.385400e+05){//debug??!!
+     fprintf (Out,"\n\n Matice K ulohy pred, time %le\n",Mp->time);
+     Smat->printmat (Out);
+     fflush(Out);//debug
+     }
+  */
   //  solution of K(r).\Delta r=\Delta F
   Mp->ssle->solve_system (Gtm,Smat,dr,fb,Out);
+  /* if (Mp->time >= 5.375400e+05 && Mp->time<= 5.385400e+05){
+     for (i=0;i<n;i++)
+     fprintf (Out,"\%le\n",dr[i]);
+     }
+  */
   //  update total displacement vector
   //  ra[j]+=dr[j];
   addv(ra, dr, n);
 ,10 → 795,14
     // re-assemble stifness matrix if it is required
     assemble_stiffness_matrix(lcid,istep,j,li,fusm);
+    //fprintf (Out,"\n\n Matice K ulohy v iteraci c. %ld, time %le\n",j,Mp->time);
+    //Smat->printmat (Out);
+    //fflush(Out);//debug
     // solution of K(r).v=F
     Mp->ssle->solve_system(Gtm,Smat,dr,fb,Out);
     // ra[k]+=dr[k];
     addv(ra, dr, n);

 /trunk/SIFEL/MEFEL/SRC/probdesc.cpp
 ,7 → 115,10
   hstrastre = (strastrestate)0;
   //  direction of unit load for homogenization for homog=7 and homog=8
   homogdir = 0;
+  //  type of macro-micro problem correspondence in homogenization
+  mami = (macromicrotype) 0;
   //  the number of Wang's tiles
   ntiletypes=0;
 ,6 → 389,13
     if (homog==9){
       fprintf (stdout,"\n general homogenization with mixed prescribed macro-strains/stresses");
     }
+    if (homog==10){
+      //  type of macro-micro problem correspondence
+      //  mami=1 - elements are connected with microproblems
+      //  mami=2 - aggregates of elements are connected with microproblems
+      xfscanf (in,"%k%ld","macmiccorres",&mami);
+      Gtm->rst=1;
+    }
   }
   // ***************************

 /trunk/SIFEL/MEFEL/SRC/seqfilesm.h
 ,6 → 45,8
 #include "global.h"
 #include "globmat.h"
 #include "graphmat.h"
+#include "homogmatm.h"
+#include "homogmech.h"
 #include "inicd.h"
 #include "intpoints.h"
 #include "j2flow.h"

 /trunk/SIFEL/MEFEL/SRC/solverm.cpp
 ,6 → 144,13
     case 7:
     case 8:{
       mechanical_homogenization();
+      // only for testing:
+      /* double *matrix_k;
+         matrix_k = new double [6];
+         paral_mechanical_homogenization (NULL, matrix_k);
+      */
       break;
     }
     default:{

 /trunk/SIFEL/METR/SRC/metrinit.cpp
 ,6 → 342,47
     //  input of loading
     Mb->read (inm);
     Tb->read (intr);
+    if ((Mp->homog==3) || (Mp->homog==5) || (Mp->homog==7) || (Mp->homog==9)){
+      //  the first order homogenization based on prescribed macro-stresses
+      //
+      // The total number of macro-stress/strain components is determined by
+      // Mt->max_ncompstr which is initialized in Mt->read function
+      // Actually, macro-sress component code numbers are supposed to be
+      // the same in all load cases => take them from the first load case
+      long *mstress_cn = Mb->give_mstress_cn(0);
+      ivector cnadd;
+      //  makerefv(cnadd, mstress_cn, Mt->max_ncompstr);
+      makerefv(cnadd, mstress_cn, Mt->max_ncompstr);
+      //  code numbers are copied from nodes to elements
+      //  and additional DOFs are added to elements
+      // do it only in the case that there are some macro-stress components prescribed (due to Mp->homog==9)
+      if (Mb->give_num_mstress_comp(0)){
+        for (long i=0; i<Mt->ne; i++)
+          Gtm->gelements[i].ndofe += Mt->give_tncomp(i);
+        Gtm->code_num_mod (cnadd, Ndofm);
+      }
+      //  volume of the domain solved
+      Mt->domvol = Mt->give_domain_vol();
+      fprintf (stdout,"\n domain volume  %le",Mt->domvol);
+    }
+    if ((Mp->homog==4) || (Mp->homog==6) || (Mp->homog==8) || (Mp->homog==9)){
+      //  the first order homogenization based on prescribed macro-strains
+      //
+      // The total number of macro-stress/strain components is determined by
+      // Mt->max_ncompstr which is initialized in Mt->read function
+      //  volume of the domain solved
+      if (Mp->homog != 9){
+        Mt->domvol = Mt->give_domain_vol();
+        fprintf (stdout,"\n domain volume  %le",Mt->domvol);
+      }
+    }
   }else{
   }

 /trunk/SIFEL/PREP/PARTITIONING/Makefile
 ,12 → 53,12
 LIBDEPS  := $(SIFEL_ROOT)$(BIN_ROOT)GEFEL/$(OUTPUTDIR)libgef.a
 LIBDEPS  += $(SIFEL_ROOT)$(BIN_ROOT)MEFEL/SRC/$(OUTPUTDIR)libmef.a
 LIBDEPS  += $(SIFEL_ROOT)$(BIN_ROOT)MEFEL/PREP/$(OUTPUTDIR)libmprep.a
-LIBDEPS  += ./metis-5.0pre2/build/Linux-i686/libmetis.a
+LIBDEPS  += ./metis-5.0pre2/build/Linux-x86_64/libmetis.a
 # path to directories with needed libraries for the preprocessors target
 DIR_LIBS := -L$(SIFEL_ROOT)$(BIN_ROOT)GEFEL/$(OUTPUTDIR)
 DIR_LIBS += -L$(SIFEL_ROOT)$(BIN_ROOT)MEFEL/SRC/$(OUTPUTDIR)
 DIR_LIBS += -L$(SIFEL_ROOT)$(BIN_ROOT)MEFEL/PREP/$(OUTPUTDIR)
-DIR_LIBS += -L./metis-5.0pre2/build/Linux-i686/
+DIR_LIBS += -L./metis-5.0pre2/build/Linux-x86_64/
 ,9 → 68,6
 all:  mechprep metispart $(MESHDECOMP)
-metispart:
-        @(cd metis-5.0pre2/; make all)
 $(MESHDECOMP): $(OBJS) $(LIBDEPS)
         @(echo "##### Creating meshdecomp executable . . .")
         $(CC) -o $@ $(OBJS) $(EXECFLAGS) $(DIR_LIBS) $(LOC_LIBS) $(SYS_LIBS)
 ,18 → 78,18
 $(OUTPUTPATH):
         @(mkdir -p $(OUTPUTPATH))
-$(LIBDEPS): mechprep
+$(LIBDEPS): mechprep metispart
 clean:
         $(RM) $(MESHDECOMP) $(OBJS) $(DEPS) *~ core
         @(echo " meshdecomp was successfully cleaned")
-cleanall: $(LIBDEPS)
+cleanall: mechprep
         @($(RM) -r -f $(addprefix $(SIFEL_ROOT)$(BIN_ROOT)$(MODULE_PATH),$(REMOVEDIR)))
         @(cd metis-5.0pre2; make realclean)
         @(echo " meshdecomp was successfully cleaned")
-cleandepall: $(LIBDEPS)
+cleandepall: mechprep
         @($(RM) -r -f $(SIFEL_ROOT)$(BIN_ROOT)$(MODULE_PATH))
         @(cd metis-5.0pre2; make realclean)
         @(echo " Preprocessors were successfully cleaned")
 ,6 → 102,8
 mechprep:
         @(cd ../../MEFEL/PREP; make $(MAKECMDGOALS))
+metispart:
+        @(cd metis-5.0pre2/; make all)
 ####################################################################
 ,4 → 119,4
 $(OUTPUTPATH)%.o : %.cpp
         @export CPATH=$(INCLUDES); \
         echo $(CC) -c $(CFLAGS) $(DEFS) $(INCLMETIS) $(<D)/$(<F); \
-        $(CC) -c $(CFLAGS) $(COMPATFLAGS) $(DEFS) $(INCLMETIS) $(<D)/$(<F) -o $(@D)/$(@F)
+        $(CC) -c $(CFLAGS) $(COMPATFLAGS) $(DEFS) $(<D)/$(<F) -o $(@D)/$(@F)

 /trunk/SIFEL/PREP/PARTITIONING/meshdecomp.cpp
 ,7 → 593,7
   FILE *out;
   long i,j;
   long adr1,adr2;
-  char outt[1025];
+  char outt[1050];
   sprintf(outt,"%s",graphf);
   out = fopen(outt,"w");
 ,7 → 680,7
   FILE *out;
   long i,j;
   long adr1,adr2;
-  char outt[1025];
+  char outt[1050];
   sprintf(outt,"%s.graph",graphf);
   out = fopen(outt,"w");
 ,7 → 761,7
   FILE *out;
   long i,j;
   long adr1,adr2;
-  char outt[1025];
+  char outt[1050];
   sprintf(outt,"%s.grf",graphf);
   out = fopen(outt,"w");
 ,7 → 832,7
   FILE *out;
   long i,j;
   long adr1,adr2;
-  char outt[1025];
+  char outt[1050];
   sprintf(outt,"%s.grf",graphf);
   out = fopen(outt,"w");
 ,7 → 962,7
   long nbn;
   long *lned,*nod,*multip,*nnsd;
   //  long **part_nodes;
-  char outt[1025];
+  char outt[2*BUFSIZ+50];
   FILE *out;
   printf("--------------------------\n");
 ,6 → 1084,9
     break;
   }
+  default:
+    fprintf(stderr, "Unknown type of mesh description (md=%d)\n", md);
+    abort();
   }
   printf ("Number of boundary nodes  %ld\n",nbn);

 /trunk/SIFEL/PREP/PARTITIONING/metis-5.0pre2/Makefile.in
 ,8 → 53,8
   #-------------------------------------------------------------------
   #Compiler information
   CC = cc
-  OPTFLAGS = -O3
-  COPTIONS = -DUNIX -pedantic -std=c99
+  OPTFLAGS = -O3
+  COPTIONS = -DUNIX -pedantic -std=c99
   #Linker information
   LDOPTIONS =

 /trunk/SIFEL/TRFEL/SRC/aliast.h
 ,8 → 118,8
 //  aliases for boundary conditions
-enum bocontypet {no_bc=0, presc_comp=1, presc_flux=2, det_climcond=3, gen_climcond=5, presc_trmiss=30, presc_trmiss_spec=31,presc_trmiss_isoheat=32,presc_trmiss_radiation=90, presc_trmiss_free_surface=40,};
-const enumstr bocontypetstr[] = {{"no_bc",0}, {"presc_comp",1}, {"presc_flux",2}, {"det_climcond",3}, {"gen_climcond",5}, {"presc_trmiss",30}, {"presc_trmiss_spec",31},{"presc_trmiss_isoheat",32},{"presc_trmiss_radiation",90},{"presc_trmiss_free_surface",40}};
+enum bocontypet {no_bc=0, presc_comp=1, presc_flux=2, det_climcond=3, gen_climcond=5, presc_trmiss=30, presc_trmiss_spec=31,presc_trmiss_isoheat=32,presc_trmiss_computed=50,presc_trmiss_radiation=90, presc_trmiss_free_surface=40,};
+const enumstr bocontypetstr[] = {{"no_bc",0}, {"presc_comp",1}, {"presc_flux",2}, {"det_climcond",3}, {"gen_climcond",5}, {"presc_trmiss",30}, {"presc_trmiss_spec",31},{"presc_trmiss_isoheat",32},{"presc_trmiss_computed",50},{"presc_trmiss_radiation",90},{"presc_trmiss_free_surface",40}};
 const kwdset bocontype_kwdset(sizeof(bocontypetstr)/sizeof(*bocontypetstr), bocontypetstr);

 /trunk/SIFEL/TRFEL/SRC/consol_wf1.cpp
 ,7 → 1144,11
     other = get_sw(pw,ipp);
     break;
   }
-  case 2:{//specific moisture content
+  case 2:{//derivative of retention curve
+    other = get_dsw_dpw(pw,ipp);
+    break;
+  }
+  case 3:{//specific moisture content
     other = get_w(pw,ipp);
     break;
   }
 ,7 → 1179,11
     fprintf (out,"Degree of saturation ()       ");
     break;
   }
-  case 2:{//moisture content
+  case 2:{//derivative of retention curve
+    fprintf (out,"Derivative of retention curve (Pa^-1)      ");
+    break;
+  }
+  case 3:{//moisture content
     fprintf (out,"Moisture (water) content (kg/kg)      ");
     //fprintf (out,"Moisture content (kg/kg)      ");
     break;

 /trunk/SIFEL/TRFEL/SRC/elemswitcht.cpp
 ,10 → 639,10
     Lqat->intpointgrad (eid);
     break;
   }
-  case ifacequadel:{
+    //case ifacequadel:{
     //Ifcquadt->intpointgrad (eid);
-    break;
-  }
+    //break;
+    //}
   case lineartett:{
     Ltett->intpointgrad (eid);
     break;
 ,24 → 685,25
     Lbt3d->intpointflux (eid);
     break;
   }
-  case barlintax:{
-    //Lbat->intpointflux (eid);
+    /* case barlintax:{
+    Lbat->intpointflux (eid);
     break;
-  }
-  case barquadt:{
-    //Qbt->intpointflux (eid);
+    }
+    case barquadt:{
+    Qbt->intpointflux (eid);
     break;
-  }
-  case barquadtax:{
+    }
+    case barquadtax:{
     //Qbat->intpointflux (eid);
     break;
-  }
+    }
+    */
   case trlint:{
-    //Ltt->intpointflux (eid);
+    Ltt->intpointflux (eid);
     break;
   }
   case trlaxisym:{
-    //Ltat->intpointflux (eid);
+    Ltat->intpointflux (eid);
     break;
   }
   case quadlint:{
 ,19 → 715,20
     break;
   }
   case quadquadtax:{
-    //Qqat->intpointflux (eid);
+    Qqat->intpointflux (eid);
     break;
   }
   case quadlaxisym:{
-    //Lqat->intpointflux (eid);
+    Lqat->intpointflux (eid);
     break;
   }
-  case ifacequadel:{
+    //     case ifacequadel:{
     //Ifcquadt->intpointflux (eid);
-    break;
-  }
+    // break;
+    // }
   case lineartett:{
-    //Ltett->intpointflux (eid);
+    Ltett->intpointflux (eid);
     break;
   }
   case linearhext:{
 ,15 → 735,14
     Lht->intpointflux (eid);
     break;
   }
-  case quadratichext:{
+    //case quadratichext:{
     //Qht->intpointflux (eid);
-    break;
-  }
+    //break;
+    //}
   default:{
     print_err("unknown element type is required",__FILE__,__LINE__,__func__);
   }
   }
 }
 /**
 ,6 → 3598,10
     Ltett->surface_flux (lcid,eid,beid,fluxes);
     break;
   }
+  case linearhext:{
+    Lht->surface_flux (lcid,eid,beid,fluxes);
+    break;
+  }
   default:{
     print_err("unknown element type is required", __FILE__, __LINE__, __func__);
   }

 /trunk/SIFEL/TRFEL/SRC/globmatt.cpp
 ,6 → 1930,13
   }
   // vector av is added into rhs vector computed in previous trfel_right_hand_side function
+  //fprintf(Outt, "%e",Tp->time);
+  //for (i=0;i<n;i++){
+  //fprintf(Outt, " %e",av[i]);
+  //}
+  //fprintf(Outt, "\n");
+  //fflush(Outt);
   addv (rhs,av,n);
   delete [] av;

 /trunk/SIFEL/TRFEL/SRC/homogmat.cpp
 ,3 → 898,18
     }
   }
 }
+/**
+   The funtion marks required non-transport quantities in the array antq.
+   @param antq - array with flags for used material types
+                 antq[i] = 1 => quantity type nontransquant(i+1) is required
+                 antq[i] = 0 => quantity type nontransquant(i+1) is not required
+   @return The function does not return anything, but it may change content of antq array.
+   TKr 01/08/2022 according to TKo
+*/
+void homogmat::give_reqntq(long *antq)
+{
+}

 /trunk/SIFEL/TRFEL/SRC/homogmat.h
 ,6 → 46,9
   void give_dof_names(namevart *dofname, long ntm);
+   /// fills array with non-transport quantities required by the model
+  void give_reqntq(long *antq);
   matrix dd;
   matrix cc;

 /trunk/SIFEL/TRFEL/SRC/homogtrans.cpp
 ,20 → 666,21
   }
   // this printing below is only for debug
-  fprintf(Outt,"\n\n vysledna matice D^M\n");
-  for (i=0;i<ntm*ncomp;i++){
-    for (j=0;j<ntm*ncomp;j++)
-      fprintf(Outt,"%e  \n",kM[i][j]);
-    fprintf(Outt,"\n");
-  }
-  fprintf(Outt,"\n\n vysledna matice C^M\n");
-  for (i=0;i<ntm;i++){
-    for (j=0;j<ntm;j++)
-      fprintf(Outt,"%e  \n",cM[i][j]);
-    fprintf(Outt,"\n");
-  }
-  fflush(Outt);
+  /* fprintf(Outt,"\n\n vysledna matice D^M\n");
+     for (i=0;i<ntm*ncomp;i++){
+     for (j=0;j<ntm*ncomp;j++)
+     fprintf(Outt,"%e  \n",kM[i][j]);
+     fprintf(Outt,"\n");
+     }
+     fprintf(Outt,"\n\n vysledna matice C^M\n");
+     for (i=0;i<ntm;i++){
+     for (j=0;j<ntm;j++)
+     fprintf(Outt,"%e  \n",cM[i][j]);
+     fprintf(Outt,"\n");
+     }
+     fflush(Outt);
+  */
   if (Kmat != NULL){
     //if (Mesprt != 0)

 /trunk/SIFEL/TRFEL/SRC/isotrmat.cpp
 ,6 → 201,10
     new_trc=0.0;
     break;
   }
+  case 50:{//heat transmission - boundary condition computed from the actual boundary flux
+    new_trc=1.0;
+    break;
+  }
   case 90:{//radiation
     new_trc=1.0;
     break;
 ,6 → 252,10
     new_trc=0.0;
     break;
   }
+  case 50:{//heat transmission - boundary condition computed from the actual boundary flux
+    new_trc=1.0;
+    break;
+  }
   case 90:{//radiation
     if(flag == 1)//for rhs vector
       new_trc=1.0;
 ,7 → 290,7
 double isotrmat::transmission_nodval(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
 {
   long k;
-  double new_nodval=0.0,t;
+  double dt=0.0,new_nodval=0.0,t;
   k=Gtt->give_dof(nn,0);
   if (k>0)   {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
 ,6 → 314,11
     new_nodval = (nodval - t);
     break;
   }
+  case 50:{//heat transmission - boundary condition computed from the actual boundary flux
+    dt = Tp->timecont.actualbacktimeincr();
+    new_nodval=Tm->ip[ipp].av[0] + Tb->fluxes[0]*nodval*dt;//cumulated value from the first integration point of the element
+    break;
+  }
   case 90:{//radiation
     //new_nodval = (nodval - t) + trc2*(nodval*nodval*nodval*nodval - t*t*t*t);//inverse eq. (trc2 = e_sigma0/trc)
     new_nodval = trc2*(nodval*nodval*nodval*nodval - t*t*t*t);//inverse eq. (trc2 = e_sigma0/trc)
 ,7 → 350,7
 double isotrmat::transmission_flux(double nodval,double trc2,long ri,long ci,long nn,long bc,long ipp)
 {
   long k;
-  double flux=0.0,t;
+  double dt=0.0,flux=0.0,t;
   k=Gtt->give_dof(nn,0);
   if (k>0)   {t = Lsrst->lhs[k-1]+Lsrst->lhsi[k-1];}
 ,6 → 366,11
     flux = (nodval - t);
     break;
   }
+  case 50:{//heat transmission - boundary condition computed from the actual boundary flux
+    dt = Tp->timecont.actualbacktimeincr();
+    flux = (Tm->ip[ipp].av[0] + Tb->fluxes[0]*nodval*dt - t);//cumulated value from the first integration point of the element
+    break;
+  }
   case 90:{//radiation
     //flux = (nodval - t) + trc2*(nodval*nodval*nodval*nodval - t*t*t*t);//inverse eq. (trc2 = e_sigma0/trc)
     flux = trc2*(nodval*nodval*nodval*nodval - t*t*t*t);//inverse eq. (trc2 = e_sigma0/trc)

 /trunk/SIFEL/TRFEL/SRC/linhext.cpp
 ,49 → 502,88
    @param eid - element id
-   JK, 12.8.2014
+   JK, 12.8.2014, revised by TKr 29/09/2022
 */
 void linhext::intpointflux (long eid)
 {
-  long i,j,k,l,ll,ii,jj,ipp;
-  vector grad(ncomp),fl(ncomp),cfl(ncomp);
-  matrix d(ncomp,ncomp);
+  long i,j,k,l,ii,jj,ipp;
-  //  loop over the number of media in the problem
   for (l=0;l<Tp->ntm;l++){
-    //  loop over all integration points on element
     for (ii=0;ii<Tp->ntm;ii++){
       for (jj=0;jj<Tp->ntm;jj++){
-        //  id of the actual integration point
-        ipp=Tt->elements[eid].ipp[ii][jj];
-        for (i=0;i<intordkm[ii][jj];i++){
-          for (j=0;j<intordkm[ii][jj];j++){
+        ipp=Tt->elements[eid].ipp[ii][jj];
+        for (i=0;i<intordkm[ii][jj];i++){
+          for (j=0;j<intordkm[ii][jj];j++){
             for(k=0;k<intordkm[ii][jj];k++){
+              //  computation of correct fluxes
+              if (Tp->fluxcomp==1)
+                Tm->computenlfluxes (l,ipp);
-              //  loop over the number of media
-              //  flux in every point is the sum of contribution from each medium
-              for (ll=0;ll<Tp->ntm;ll++){
-                //  gradients
-                Tm->givegrad (ll,ipp,grad);
-                //  block of conductivity matrix of material
-                Tm->matcond (d,ipp,l,ll);
-                //  flux caused by the actual gradient
-                mxv (d,grad,fl);
-                //  summation of flux contributions
-                addv (fl,cfl,cfl);
-              }
-              Tm->storeflux (l,ipp,cfl);
-              ipp++;
-            }
-          }
-        }
+              ipp++;
+            }
+          }
+        }
+        for (i=0;i<intordcm[ii][jj];i++){
+          for (j=0;j<intordcm[ii][jj];j++){
+            for(k=0;k<intordcm[ii][jj];k++){
+              //  computation of correct fluxes
+              if (Tp->fluxcomp==1)
+                Tm->computenlfluxes (l,ipp);
+              ipp++;
+            }
+          }
+        }
+        if (Tp->savemode==1)  break;
       }
-      if (Tp->savemode==1)
-        break;
+      if (Tp->savemode==1)  break;
     }
   }
+  /*
+      long i,j,k,l,ll,ii,jj,ipp;
+      vector grad(ncomp),fl(ncomp),cfl(ncomp);
+      matrix d(ncomp,ncomp);
+      //  loop over the number of media in the problem
+      for (l=0;l<Tp->ntm;l++){
+      //  loop over all integration points on element
+      for (ii=0;ii<Tp->ntm;ii++){
+      for (jj=0;jj<Tp->ntm;jj++){
+      //  id of the actual integration point
+      ipp=Tt->elements[eid].ipp[ii][jj];
+      for (i=0;i<intordkm[ii][jj];i++){
+      for (j=0;j<intordkm[ii][jj];j++){
+      for(k=0;k<intordkm[ii][jj];k++){
+      //  loop over the number of media
+      //  flux in every point is the sum of contribution from each medium
+      for (ll=0;ll<Tp->ntm;ll++){
+      //  gradients
+      Tm->givegrad (ll,ipp,grad);
+      //  block of conductivity matrix of material
+      Tm->matcond (d,ipp,l,ll);
+      //  flux caused by the actual gradient
+      mxv (d,grad,fl);
+      //  summation of flux contributions
+      addv (fl,cfl,cfl);
+      }
+      Tm->storeflux (l,ipp,cfl);
+      ipp++;
+      }
+      }
+      }
+      }
+      if (Tp->savemode==1)
+      break;
+      }
+      }
+  */
 }
 /**
 ,3 → 2664,65
   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 26/09/2022 according to JK
+*/
+void linhext::surface_flux (long lcid,long eid,long beid,double *fluxes)
+{
+  long i,ipp,fid;
+  double q,qn,area;
+  bocontypet *bc;
+  vector x(ASTCKVEC(nne)),y(ASTCKVEC(nne)),z(ASTCKVEC(nne));
+  vector bb(ASTCKVEC(4)),cc(ASTCKVEC(4)),dd(ASTCKVEC(4)),volcoord(ASTCKVEC(4));
+  vector r(ASTCKVEC(ndofe)),t(ASTCKVEC(nne)),grad(ASTCKVEC(ncomp)),flux(ASTCKVEC(ncomp)),fl(ASTCKVEC(ncomp)),n(ASTCKVEC(ncomp));
+  matrix gm(ASTCKMAT(ncomp,nne)),d(ASTCKMAT(ncomp,ncomp));
+  //  node coordinates
+  Tt->give_node_coord3d (x,y,z,eid);
+  if (Tp->savemode==0)
+    ipp=Tt->elements[eid].ipp[lcid][lcid];
+  if (Tp->savemode==1)
+    ipp=Tt->elements[eid].ipp[0][0];
+  //the first int. point
+  Tm->computenlfluxes (lcid,ipp);
+  Tm->givefluxes (lcid,ipp,flux);
+  //  indicators of boundary conditions
+  bc = new bocontypet [nsurf];
+  Tb->bf[lcid].elemload[beid].give_bc (bc);
+  //  loop over surfaces
+  for (i=0;i<nsurf;i++){
+    if (bc[i]==1){
+      //  function constructs the outer normal vector
+      hexahedra_normal_vectors (i,x,y,z,n);
+      //  area of the i-th element surface
+      area = hexahedra_surface_areas (i,x,y,z);
+      //  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/linhext.h
43,6 → 43,7
void res_internal_fluxes (long eid,vector &elemif);
double total_integral (long eid,vector &nodval);
void res_boundary_flux (vector &f,long lcid,long eid,long leid);
void surface_flux (long lcid,long eid,long beid,double *fluxes);


void nod_grads_ip (long eid);

 /trunk/SIFEL/TRFEL/SRC/lintett.cpp
 ,49 → 1488,79
    @param eid - element id
-   JK, 12.8.2014
+   JK, 12.8.2014, revised by TKr 29/09/2022
 */
 void lintett::intpointflux (long eid)
 {
-  long i,j,k,l,ll,ii,jj,ipp;
-  vector grad(ASTCKVEC(ncomp)),fl(ASTCKVEC(ncomp)),cfl(ASTCKVEC(ncomp));
-  matrix d(ASTCKMAT(ncomp,ncomp));
+  long i,l,ii,jj,ipp;
-  //  loop over the number of media in the problem
   for (l=0;l<Tp->ntm;l++){
-    //  loop over all integration points on element
     for (ii=0;ii<Tp->ntm;ii++){
       for (jj=0;jj<Tp->ntm;jj++){
-        //  id of the actual integration point
-        ipp=Tt->elements[eid].ipp[ii][jj];
-        for (i=0;i<intordkm[ii][jj];i++){
-          for (j=0;j<intordkm[ii][jj];j++){
-            for(k=0;k<intordkm[ii][jj];k++){
-              //  loop over the number of media
-              //  flux in every point is the sum of contribution from each medium
-              for (ll=0;ll<Tp->ntm;ll++){
-                //  gradients
-                Tm->givegrad (ll,ipp,grad);
-                //  block of conductivity matrix of material
-                Tm->matcond (d,ipp,l,ll);
-                //  flux caused by the actual gradient
-                mxv (d,grad,fl);
-                //  summation of flux contributions
-                addv (fl,cfl,cfl);
-              }
-              Tm->storeflux (l,ipp,cfl);
-              ipp++;
-            }
-          }
-        }
+        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 (l,ipp);
+          ipp++;
+        }
+        for (i=0;i<intordcm[ii][jj];i++){
+          //  computation of correct fluxes
+          if (Tp->fluxcomp==1)
+            Tm->computenlfluxes (l,ipp);
+          ipp++;
+        }
+        if (Tp->savemode==1)  break;
       }
-      if (Tp->savemode==1)
-        break;
+      if (Tp->savemode==1)  break;
     }
   }
+  /*
+      long i,j,k,l,ll,ii,jj,ipp;
+      vector grad(ASTCKVEC(ncomp)),fl(ASTCKVEC(ncomp)),cfl(ASTCKVEC(ncomp));
+      matrix d(ASTCKMAT(ncomp,ncomp));
+      //  loop over the number of media in the problem
+      for (l=0;l<Tp->ntm;l++){
+      //  loop over all integration points on element
+      for (ii=0;ii<Tp->ntm;ii++){
+      for (jj=0;jj<Tp->ntm;jj++){
+      //  id of the actual integration point
+      ipp=Tt->elements[eid].ipp[ii][jj];
+      for (i=0;i<intordkm[ii][jj];i++){
+      for (j=0;j<intordkm[ii][jj];j++){
+      for(k=0;k<intordkm[ii][jj];k++){
+      //  loop over the number of media
+      //  flux in every point is the sum of contribution from each medium
+      for (ll=0;ll<Tp->ntm;ll++){
+      //  gradients
+      Tm->givegrad (ll,ipp,grad);
+      //  block of conductivity matrix of material
+      Tm->matcond (d,ipp,l,ll);
+      //  flux caused by the actual gradient
+      mxv (d,grad,fl);
+      //  summation of flux contributions
+      addv (fl,cfl,cfl);
+      }
+      Tm->storeflux (l,ipp,cfl);
+      ipp++;
+      }
+      }
+      }
+      }
+      if (Tp->savemode==1)
+      break;
+      }
+      }
+  */
 }
 ,7 → 2490,7
   //  id of integration point for determination of the D matrix
   ipp=Tt->elements[eid].ipp[0][0];
-  //  determinant fr the volume coordinates
+  //  determinant for the volume coordinates
   det = det3d (x.a,y.a,z.a);
   if(det < 0.0){

 /trunk/SIFEL/TRFEL/SRC/npsolvert.cpp
 ,7 → 566,6
   //Kmat->copygm (*Cmat);
   //new:
   Cmat->addgm (dt*alpha,*Kmat);
   //fprintf(Outt, "\nOverall matrix at istep=%ld, jstep=%ld:\n", Tp->istep, Tp->jstep);
   //Cmat->printmat(Outt);
 ,9 → 573,11
   trfel_right_hand_side (0,rhs,n);
   //fprintf(Outt, "\nRight-hand side vector:\n");
+  //  fprintf(Outt, "%e",Tp->time);
   //for (j=0;j<n;j++){
-  //fprintf(Outt, "%e\n",rhs[j]);
+  //fprintf(Outt, " %e",rhs[j]);
   //}
+  //fprintf(Outt, "\n");
   //fflush(Outt);
   //mechanical influence of transport problem
 ,9 → 584,11
   trfel_right_hand_side2 (0,rhs,n);
   //fprintf(Outt, "\nRight-hand2 side vector:\n");
+  //fprintf(Outt, "%e",Tp->time);
   //for (j=0;j<n;j++){
-  //fprintf(Outt, "%e\n",rhs[j]);
+  //fprintf(Outt, " %e",rhs[j]);
   //}
+  //fprintf(Outt, "\n");
   //fflush(Outt);
   //  computation of the right hand side vector

 /trunk/SIFEL/TRFEL/SRC/quadlinaxisym.cpp
 ,7 → 1310,7
 */
 void quadlinaxisym::intpointflux (long eid)
 {
-  long i,ii,jj,k,ipp;
+  long i,j,ii,jj,k,ipp;
   for (k=0;k<Tp->ntm;k++){
 ,20 → 1319,23
         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<intordkm[ii][jj];i++){
+          for (j=0;j<intordkm[ii][jj];j++){
+            //  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++;
+        for (i=0;i<intordcm[ii][jj];i++){
+          for (j=0;j<intordcm[ii][jj];j++){
+            //  computation of correct fluxes
+            if (Tp->fluxcomp==1)
+              Tm->computenlfluxes (k,ipp);
+            ipp++;
+          }
         }
         if (Tp->savemode==1)  break;
       }

 /trunk/SIFEL/TRFEL/SRC/quadlineart.cpp
 ,7 → 2657,7
 */
 void quadlineart::intpointflux (long eid)
 {
-  long i,ii,jj,k,ipp;
+  long i,j,ii,jj,k,ipp;
   for (k=0;k<Tp->ntm;k++){
 ,20 → 2666,23
         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<intordkm[ii][jj];i++){
+          for (j=0;j<intordkm[ii][jj];j++){
+            //  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++;
+        for (i=0;i<intordcm[ii][jj];i++){
+          for (j=0;j<intordcm[ii][jj];j++){
+            //  computation of correct fluxes
+            if (Tp->fluxcomp==1)
+              Tm->computenlfluxes (k,ipp);
+            ipp++;
+          }
         }
         if (Tp->savemode==1)  break;
       }
 ,7 → 2826,49
 }
+/**
+   function computes other values in nodes of element
+   @param eid - element id
+   TKr, 05/05/2022 according to TKo nod_eqother_ip (long eid)
+*/
+void quadlineart::nod_other_ip (long eid)
+{
+  long i,j,k,ipp,ncompo;
+  double area;
+  ivector ipnum(ASTCKIVEC(nne)), nod(ASTCKIVEC(nne));
+  vector other;
+  //  numbers of integration points closest to nodes
+  //  (function is from the file GEFEL/ordering.cpp)
+  ipp=Tt->elements[eid].ipp[0][0];
+  nodip_planelq (ipp,intordkm[0][0],ipnum);
+  //  node numbers of the element
+  Tt->give_elemnodes (eid,nod);
+  for (i=0;i<nne;i++){
+    ncompo = Tm->givencompother ();
+    reallocv (RSTCKVEC(ncompo,other));
+    Tm->giveother (ipnum[i],0,ncompo,other.a);
+    //  storage of other components to the node
+    j=nod[i];
+    for (k=0; k<ncompo; k++){
+      if (Tp->otheraver==0 || Tp->otheraver==1)
+        Tt->nodes[j].storeother (k,other(k));
+      if (Tp->otheraver==2){
+        area = element_area (eid);
+        Tt->nodes[j].storeother (k,area,other(k));
+      }
+    }
+  }
+}
 /**
    function computes other values directly in nodes of element

/trunk/SIFEL/TRFEL/SRC/quadlineart.h
74,6 → 74,7
void nod_grads_ip (long eid);
void nod_fluxes_ip (long eid);
void nod_eqother_ip (long eid);
void nod_other_ip (long eid);
void nod_others_comp (long lcid,long eid,long ri,long ci);

void higher_to_lower_level (long eid,long counter,double buff);

 /trunk/SIFEL/TRFEL/SRC/quadquadrilatt.cpp
 ,7 → 1557,7
 */
 void quadquadrilatt::intpointflux (long eid)
 {
-  long i,ii,jj,k,ipp;
+  long i,j,ii,jj,k,ipp;
   for (k=0;k<Tp->ntm;k++){
 ,20 → 1566,23
         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<intordkm[ii][jj];i++){
+          for (j=0;j<intordkm[ii][jj];j++){
+            //  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++;
+        for (i=0;i<intordcm[ii][jj];i++){
+          for (j=0;j<intordcm[ii][jj];j++){
+            //  computation of correct fluxes
+            if (Tp->fluxcomp==1)
+              Tm->computenlfluxes (k,ipp);
+            ipp++;
+          }
         }
         if (Tp->savemode==1)  break;
       }

 /trunk/SIFEL/TRFEL/SRC/quadquadrilattax.cpp
 ,7 → 1190,52
+/**
+   function computes correct fluxes at integration points on element
+   @param eid - element id
+   TKr, 29/09/2022
+*/
+void quadquadrilattax::intpointflux (long eid)
+{
+  long i,j,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++){
+          for (j=0;j<intordkm[ii][jj];j++){
+            //  computation of correct fluxes
+            if (Tp->fluxcomp==1)
+              Tm->computenlfluxes (k,ipp);
+            ipp++;
+          }
+        }
+        for (i=0;i<intordcm[ii][jj];i++){
+          for (j=0;j<intordcm[ii][jj];j++){
+            //  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 other values in nodes of element

 /trunk/SIFEL/TRFEL/SRC/quadquadrilattax.h
 ,6 → 41,8
   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 intpointflux (long eid);
   void nod_others (long lcid,long eid,long ri,long ci);
   void convection_vector (vector &v,long lcid,long eid,long leid,long ri,long ci);

 /trunk/SIFEL/TRFEL/SRC/transmat.cpp
 ,7 → 2470,7
   switch (Tp->mednam)//names of transported media
     {
     case water:{
-      ncompo=3;
+      ncompo=4;
       break;
     }
     case moisture:{
 ,9 → 3694,12
         moisth[i].give_reqntq(antq);
         break;
   }
+  case homomat:{
+    hommat[i].give_reqntq(antq);
+    break;
+  }
     /* case nlisotransmat:
-       case homomat:
        case discontisotrmat:
        case cernyconcrete:
        case bazantpedersen:

 /trunk/SIFEL/TRFEL/SRC/trlinaxisym.cpp
 ,7 → 1218,50
 /**
+   function computes correct fluxes at integration points on element
+   @param eid - element id
+   TKr, 01/02/2010
+*/
+void trlinaxisym::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 other values in nodes of element
    @param lcid - load case id

/trunk/SIFEL/TRFEL/SRC/trlinaxisym.h
45,6 → 45,7
void res_volume_rhs_vector (vector &f,long eid,long lcid);
void res_volume_rhs_vector2 (vector &f,long eid,long lcid);

void intpointflux (long eid);

void nod_others (long lcid,long eid,long ri,long ci);
void nod_others_comp (long lcid,long eid,long ri,long ci);

 /trunk/SIFEL/TRFEL/SRC/van_genuchten.cpp
 ,6 → 11,7
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
+#include "globalt.h"
 #include "van_genuchten.h"
 van_genuchten_reten::van_genuchten_reten()
 ,7 → 19,7
   ssat = 1.0;      //saturation degree [-]
   sirr = 0.2;      //irreversible saturation degree [-]
   ksat = 0.1;      //saturation permeability [??]
-  gamaw = 10.0; //water weigth [N/m^3]
+  gamaw = 10.0;    //water weigth [kN/m^3]
   delta = 20.0;    //parameter in m^(-1); delta = 20.0 [1/m]; delta = 0.2 [1/cm]
   expn = 2.0;      //parameter n
   expm = 0.0;      //parameter m
 ,8 → 65,8
 {
   double sw,hp,theta;
-  pw = pw/1000.0;
+  pw = pw/1000.0;//in kPa
   if(pw >= 0.0)
     pw = 0.0;//for positive pressure
   hp = pw/gamaw;
 ,10 → 89,10
 */
 double van_genuchten_reten::dsw_dpw(double pw)
 {
-  double dsw_dpw,hp;
+  double dsw_dpw,hp,expm;
-  pw = pw/1000.0;
+  pw = pw/1000.0;//in kPa
   if(pw >= 0.0){
     pw = 0.0;//for positive pressure
     dsw_dpw = 0.0;
 ,7 → 103,7
       expm = 1.0 - 1.0/expn;
     dsw_dpw = (expm*expn*delta*(ssat - sirr)*pow(delta*fabs(hp),(expn-1)))/(pow(1+(pow(delta*fabs(hp),expn)),(expm+1)))/gamaw;
   }
   dsw_dpw =  dsw_dpw/1000.0;//basic units are Pa
   return(dsw_dpw);

Subversion Repositories sifel

Compare Revisions

Ignore whitespace Rev 775 → Rev 776