Starting implementing mixture laws

include/micro.hpp

@@ -158,12 +158,8 @@ enum {
 enum {
        	LINEAR,
        	ONE_WAY,
-       	FULL
-};
-
-enum {
-	HOMOG_LINEAR,
-	HOMOG_NON_LINEAR
+       	FULL,
+       	RULE_MIXTURE_1
 };
 
 
@@ -184,7 +180,7 @@ class micropp {
 		const double lx, ly, lz;
 		const double dx, dy, dz;
 		const double vol_tot;
-		const double wg, ivol;
+		const double wg, ivol, evol;
 
 		const int micro_type, nvars;
 		const int nsubiterations;
@@ -230,6 +226,10 @@ class micropp {
 		int its_with_A0;
 		ell_matrix *A0;
 
+		/* Rule of Mixture Stuff (for 2 mats micro-structure only) */
+		double Vm;  // Volume fraction of Matrix
+		double Vf;  // Volume fraction of Fiber
+
 		/* IO files */
 		const bool write_log_flag;
 		int log_id = 0;
@@ -243,6 +243,7 @@ class micropp {
 
 		void homogenize_linear(gp_t<tdim> *gp_ptr);
 		void homogenize_non_linear(gp_t<tdim> *gp_ptr);
+		void homogenize_rule_mixture_1(gp_t<tdim> *gp_ptr);
 
 		void calc_ctan_lin();
 
@@ -282,6 +283,8 @@ class micropp {
 
 		void calc_bmat(int gp, double bmat[nvoi][npe * dim]) const;
 
+		void calc_volume_fractions();
+
 		bool calc_vars_new(const double *u, const double *vars_old,
 				   double *vars_new) const;
 
@@ -328,7 +331,7 @@ class micropp {
 
 		void get_ctan(const int gp_id, double *ctan) const;
 
-		void homogenize(const int homog_type = HOMOG_NON_LINEAR);
+		void homogenize();
 
 		/* Extras */
 

include/micropp_c.h

@@ -49,7 +49,7 @@ extern "C" {
 	void micropp3_get_ctan(const struct micropp3 *self, const int gp_id,
 			       double *ctan);
 
-	void micropp3_homogenize(struct micropp3 *self, const int homog_type);
+	void micropp3_homogenize(struct micropp3 *self);
 
 	void micropp3_update_vars(struct micropp3 *self);
 

src/homogenize.cpp

@@ -56,7 +56,7 @@ void micropp<tdim>::get_ctan(const int gp_id, double *ctan) const
 
 
 template <int tdim>
-void micropp<tdim>::homogenize(const int homog_type)
+void micropp<tdim>::homogenize()
 {
 	INST_START;
 
@@ -65,8 +65,11 @@ void micropp<tdim>::homogenize(const int homog_type)
 
 		gp_t<tdim> *gp_ptr = &gp_list[igp];
 
-		if ((homog_type == HOMOG_LINEAR && gp_ptr->coupling == ONE_WAY)
-		    || gp_ptr->coupling == LINEAR) {
+		if (gp_ptr->coupling == RULE_MIXTURE_1) {
+
+			homogenize_rule_mixture_1(gp_ptr);
+
+		} else if (gp_ptr->coupling == LINEAR) {
 
 			/*
 			 * Computational cheap calculation
@@ -94,7 +97,8 @@ void micropp<tdim>::homogenize_linear(gp_t<tdim> * gp_ptr)
 	memset (gp_ptr->stress, 0.0, nvoi * sizeof(double));
 	for (int i = 0; i < nvoi; ++i) {
 		for (int j = 0; j < nvoi; ++j) {
-			gp_ptr->stress[i] += ctan_lin[i * nvoi + j] * gp_ptr->strain[j];
+			gp_ptr->stress[i] += 
+				ctan_lin[i * nvoi + j] * gp_ptr->strain[j];
 		}
 	}
 }
@@ -211,6 +215,13 @@ void micropp<tdim>::homogenize_non_linear(gp_t<tdim> * gp_ptr)
 }
 
 
+template<int tdim>
+void micropp<tdim>::homogenize_rule_mixture_1(gp_t<tdim> * gp_ptr)
+{
+	
+}
+
+
 template <int tdim>
 void micropp<tdim>::update_vars()
 {

src/micro.f95

@@ -85,12 +85,10 @@ module libmicropp
        real(c_double), intent(out), dimension(*) :: ctan
      end subroutine micropp3_get_ctan
 
-     subroutine micropp3_homogenize(this, homog_type) bind(C)
-       use, intrinsic :: iso_c_binding, only: c_int
+     subroutine micropp3_homogenize(this) bind(C)
        import micropp3
        implicit none
        type(micropp3), intent(inout) :: this
-       integer(c_int), intent(in), value :: homog_type
      end subroutine micropp3_homogenize
 
      logical(c_bool) function micropp3_is_non_linear(this, gp_id) bind(C)

src/micro_common.cpp

@@ -21,6 +21,7 @@
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  */
 
+
 #include "micro.hpp"
 #include "material.hpp"
 
@@ -49,6 +50,7 @@ micropp<tdim>::micropp(const micropp_params_t &params):
 	wg(((tdim == 3) ? dx * dy * dz : dx * dy) / npe),
 	vol_tot((tdim == 3) ? lx * ly * lz : lx * ly),
 	ivol(1.0 / (wg * npe)),
+	evol((tdim == 3) ? dx * dy * dz : dx * dy),
 	micro_type(params.type), nvars(nelem * npe * NUM_VAR_GP),
 
 	nr_max_its(params.nr_max_its),
@@ -126,6 +128,8 @@ micropp<tdim>::micropp(const micropp_params_t &params):
 		}
 	}
 
+	calc_volume_fractions();
+
 	if (params.use_A0) {
 #ifdef _OPENMP
 		int num_of_A0s = omp_get_max_threads();
@@ -815,8 +819,10 @@ void micropp<tdim>::print_info() const
 			cout << "NO TYPE" << endl;
 			break;
 	}
+	cout << "MATRIX [%]  : " << Vm << endl;
+	cout << "FIBER  [%]  : " << Vf << endl;
 
-	cout << "LINEAR : " << num_no_coupling << " GPs" << endl;
+	cout << "LINEAR      : " << num_no_coupling << " GPs" << endl;
 	cout << "ONE_WAY     : " << num_one_way     << " GPs" << endl;
 	cout << "FULL        : " << num_full        << " GPs" << endl;
 	cout << "USE A0      : " << use_A0 << endl;
@@ -907,6 +913,28 @@ void micropp<tdim>::calc_ave_stress(const double *u, double stress_ave[nvoi],
 }
 
 
+template <int tdim>
+void micropp<tdim>::calc_volume_fractions()
+{
+	Vm = 0.0;
+	Vf = 0.0;
+	for (int ez = 0; ez < nez; ++ez) {
+		for (int ey = 0; ey < ney; ++ey) {
+			for (int ex = 0; ex < nex; ++ex) {
+				const int e_i = glo_elem(ex, ey, ez);
+				if (elem_type[e_i] == 0) {
+					Vm += evol;
+				} else if (elem_type[e_i] >= 1) {
+					Vf += evol;
+				}
+			}
+		}
+	}
+	Vm /= vol_tot;
+	Vf /= vol_tot;
+}
+
+
 template <int tdim>
 void micropp<tdim>::calc_ave_strain(const double *u, double strain_ave[nvoi]) const
 {

src/micropp_c.cpp

@@ -82,10 +82,10 @@ extern "C" {
 		ptr->get_ctan(gp, ctan);
 	}
 
-	void micropp3_homogenize(micropp3 *self, const int homog_type)
+	void micropp3_homogenize(micropp3 *self)
 	{
 		micropp<3> *ptr = (micropp<3> *) self->ptr;
-		ptr->homogenize(homog_type);
+		ptr->homogenize();
 	}
 
 	int micropp3_get_cost(const micropp3 *self, int gp_id)

test/benchmark-mic-1.cpp

@@ -84,10 +84,10 @@ int main(int argc, char **argv)
 	mic_params.calc_ctan_lin = false;
 	mic_params.lin_stress = false;
 
-	mic_params.print();
+	//mic_params.print();
 
 	micropp<3> micro(mic_params);
-	//micro.print_info();
+	micro.print_info();
 
 	if (print) {
 		char filename[128];

test/test3d_3.f90

@@ -101,7 +101,7 @@ program test3d_3
         end if
 
         call micropp3_set_strain(micro, gp_id, eps)
-        call micropp3_homogenize(micro, 1)
+        call micropp3_homogenize(micro)
 
         call micropp3_get_stress(micro, gp_id, sig)
         call micropp3_get_ctan(micro, gp_id, ctan)

test/test3d_6.c

@@ -76,7 +76,7 @@ int main (int argc, char *argv[])
 			eps[dir] += D_EPS;
 
 		micropp3_set_strain(&micro, 0, eps);
-		micropp3_homogenize(&micro, 1);
+		micropp3_homogenize(&micro);
 		micropp3_get_stress(&micro, 0, sig);
 		micropp3_get_ctan(&micro, 0, ctan);
 		int sigma_cost = micropp3_get_cost(&micro, 0);