Adding benchmark and MIXTURE RULE CHAMIS

include/micro.hpp

@@ -161,7 +161,7 @@ enum {
        	FE_LINEAR,
        	FE_ONE_WAY,
        	FE_FULL,
-       	RULE_MIXTURE_LIN_1
+       	MIXTURE_RULE_CHAMIS
 };
 
 
@@ -238,24 +238,25 @@ class micropp {
 		int log_id = 0;
 		ofstream ofstream_log;
 
-		/* GPU number of device selection */
+		/* GPU number for device selection */
 		int gpu_id = 0;
 
 
 		/* Private function members */
 
-		/* Linear homogenizations */
+		/*
+		 * Linear homogenizations 
+		 * Applies to FE RVE model and Mixture rules
+		 *
+		 */
 		void homogenize_linear(gp_t<tdim> *gp_ptr);
 
 		/* FE-based homogenizations */
 		void homogenize_fe_one_way(gp_t<tdim> *gp_ptr);
 		void homogenize_fe_full(gp_t<tdim> *gp_ptr);
 
-		/* Rule of mixture (cheap) */
-		void homogenize_rule_mixture_1(gp_t<tdim> *gp_ptr);
-
 		void calc_ctan_lin_fe_models();
-		void calc_ctan_lin_rule_mixture_lin_1(double ctan[nvoi * nvoi]);
+		void calc_ctan_lin_mixture_rule_Chamis(double ctan[nvoi * nvoi]);
 
 		material_t *get_material(const int e) const;
 

src/homogenize.cpp

@@ -87,11 +87,16 @@ void micropp<tdim>::homogenize()
 
 		gp_t<tdim> *gp_ptr = &gp_list[igp];
 
-		if (gp_ptr->coupling == FE_LINEAR || gp_ptr->coupling == RULE_MIXTURE_LIN_1) {
+		if (gp_ptr->coupling == FE_LINEAR ||
+		    gp_ptr->coupling == MIXTURE_RULE_CHAMIS) {
 
 			/*
 			 * Computational cheap calculation
 			 * stress = ctan_lin * strain
+			 *
+			 * All mixture rules are linear in Micropp
+			 * so the homogenization of the stress tensor
+			 * is this simple and cheap procedure.
 			 */
 
 			homogenize_linear(gp_ptr);
@@ -317,13 +322,6 @@ void micropp<tdim>::homogenize_fe_full(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_common.cpp

@@ -103,8 +103,7 @@ micropp<tdim>::micropp(const micropp_params_t &params):
 		gp_list[gp].nvars = nvars;
 
 		if (params.coupling == nullptr ||
-		    (params.coupling[gp] == FE_ONE_WAY ||
-		     params.coupling[gp] == FE_FULL)) {
+		    (params.coupling[gp] == FE_ONE_WAY || params.coupling[gp] == FE_FULL)) {
 			gp_list[gp].allocate_u();
 		}
 	}
@@ -160,15 +159,16 @@ micropp<tdim>::micropp(const micropp_params_t &params):
 	}
 
 	for (int gp = 0; gp < ngp; ++gp) {
+
 		if (gp_list[gp].coupling == FE_LINEAR || gp_list[gp].coupling == FE_ONE_WAY ||
 		    gp_list[gp].coupling == FE_FULL) {
 
 			memcpy(gp_list[gp].ctan, ctan_lin_fe, nvoi * nvoi * sizeof(double));
 
-		} else if (gp_list[gp].coupling == RULE_MIXTURE_LIN_1) {
+		} else if (gp_list[gp].coupling == MIXTURE_RULE_CHAMIS) {
 
 			double ctan[nvoi * nvoi];
-			calc_ctan_lin_rule_mixture_lin_1(ctan);
+			calc_ctan_lin_mixture_rule_Chamis(ctan);
 			memcpy(gp_list[gp].ctan, ctan, nvoi * nvoi * sizeof(double));
 
 		}
@@ -261,9 +261,11 @@ template <int tdim>
 int micropp<tdim>::get_non_linear_gps(void) const
 {
 	int count = 0;
-	for (int gp = 0; gp < ngp; ++gp)
-		if (gp_list[gp].allocated)
+	for (int gp = 0; gp < ngp; ++gp) {
+		if (gp_list[gp].allocated) {
 			count ++;
+		}
+	}
 	return count;
 }
 
@@ -304,8 +306,63 @@ void micropp<tdim>::calc_ctan_lin_fe_models()
 
 
 template <int tdim>
-void micropp<tdim>::calc_ctan_lin_rule_mixture_lin_1(double ctan[nvoi * nvoi])
+void micropp<tdim>::calc_ctan_lin_mixture_rule_Chamis(double ctan[nvoi * nvoi])
 {
+
+ const double Em = material_list[0]->E;
+ const double nu_m = material_list[0]->nu;
+ const double Gm = Em / (2 * (1 + nu_m));
+
+ const double Ef = material_list[1]->E;
+ const double nu_f = material_list[1]->nu;
+ const double Gf = Ef / (2 * (1 + nu_f));
+
+ const double E11 = Vf * Ef + Vm * Em;
+ const double E22 = Em / (1 - sqrt(Vf) * (1 - Em / Ef));
+ const double nu_12 = Vf * nu_f + Vm * nu_m;
+ const double G12 = Gm / (1 - sqrt(Vf) * (1 - Gm / Gf));
+
+ const double S[3][3] = {
+	 {       1 / E11, - nu_12 / E11, - nu_12 / E11 },
+	 { - nu_12 / E11,       1 / E22, - nu_12 / E22 },
+	 { - nu_12 / E11, - nu_12 / E22,       1 / E22 },
+ };
+
+ double det =
+	 S[0][0] * (S[1][1] * S[2][2] - S[2][1] * S[1][2]) -
+	 S[0][1] * (S[1][0] * S[2][2] - S[2][0] * S[1][2]) +
+	 S[0][2] * (S[1][0] * S[2][1] - S[2][0] * S[1][1]);
+
+ double c00 = +(S[1][1] * S[2][2] - S[2][1] * S[1][2]);
+ double c01 = -(S[1][0] * S[2][2] - S[2][0] * S[1][2]);
+ double c02 = +(S[1][0] * S[2][1] - S[2][0] * S[1][1]);
+
+ double c10 = -(S[0][1] * S[2][2] - S[2][1] * S[0][2]);
+ double c11 = +(S[0][0] * S[2][2] - S[2][0] * S[0][2]);
+ double c12 = -(S[0][0] * S[2][1] - S[2][0] * S[0][1]);
+
+ double c20 = +(S[0][1] * S[1][2] - S[1][1] * S[0][2]);
+ double c21 = -(S[0][0] * S[1][2] - S[1][0] * S[0][2]);
+ double c22 = +(S[0][0] * S[1][1] - S[1][0] * S[0][1]);
+
+ memset (ctan, 0, nvoi * nvoi * sizeof(double));
+
+ ctan[0 * nvoi + 0] = c00 / det;
+ ctan[0 * nvoi + 1] = c01 / det;
+ ctan[0 * nvoi + 2] = c02 / det;
+
+ ctan[1 * nvoi + 0] = c10 / det;
+ ctan[1 * nvoi + 1] = c11 / det;
+ ctan[1 * nvoi + 2] = c12 / det;
+
+ ctan[2 * nvoi + 0] = c12 / det;
+ ctan[2 * nvoi + 1] = c21 / det;
+ ctan[2 * nvoi + 2] = c22 / det;
+
+ ctan[3 * nvoi + 3] = G12;
+ ctan[4 * nvoi + 4] = G12;
+ ctan[5 * nvoi + 5] = G12;
+
 }
 
 
@@ -875,12 +932,14 @@ void micropp<tdim>::calc_ave_stress(const double *u, double stress_ave[nvoi],
 					double stress_gp[nvoi], strain_gp[nvoi];
 					get_strain(u, gp, strain_gp, ex, ey, ez);
 					get_stress(gp, strain_gp, vars_old, stress_gp, ex, ey, ez);
-					for (int v = 0; v < nvoi; ++v)
+					for (int v = 0; v < nvoi; ++v) {
 						stress_aux[v] += stress_gp[v] * wg;
+					}
 
 				}
-				for (int v = 0; v < nvoi; ++v)
+				for (int v = 0; v < nvoi; ++v) {
 					stress_ave[v] += stress_aux[v];
+				}
 			}
 		}
 	}
@@ -927,18 +986,21 @@ void micropp<tdim>::calc_ave_strain(const double *u, double strain_ave[nvoi]) co
 					double strain_gp[nvoi];
 
 					get_strain(u, gp, strain_gp, ex, ey, ez);
-					for (int v = 0; v < nvoi; ++v)
+					for (int v = 0; v < nvoi; ++v) {
 						strain_aux[v] += strain_gp[v] * wg;
+					}
 				}
 
-				for (int v = 0; v < nvoi; v++)
+				for (int v = 0; v < nvoi; v++) {
 					strain_ave[v] += strain_aux[v];
+				}
 			}
 		}
 	}
 
-	for (int v = 0; v < nvoi; v++)
+	for (int v = 0; v < nvoi; v++) {
 		strain_ave[v] /= vol_tot;
+	}
 }
 
 

src/micropp_c.cpp

@@ -53,6 +53,8 @@ extern "C" {
 		params.write_log = true;
 
 		self->ptr = new micropp<3>(params);
+
+		delete [] params.coupling;
 	}
 
 	void micropp3_free(micropp3 *self)

test/CMakeLists.txt

@@ -25,6 +25,7 @@ set(testsources
 	benchmark-mic-1.cpp
 	benchmark-mic-2.cpp
 	benchmark-mic-3.cpp
+	benchmark-mic-4.cpp
 	)
 
 # Iterate over the list above

test/benchmark-mic-4.cpp

+/*
+ *  This is a test example for MicroPP: a finite element library
+ *  to solve microstructural problems for composite materials.
+ *
+ *  Copyright (C) - 2018 - Jimmy Aguilar Mena <kratsbinovish@gmail.com>
+ *                         Guido Giuntoli <gagiuntoli@gmail.com>
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+
+#include <iostream>
+#include <iomanip>
+#include <cstring>
+#include <cassert>
+#include <chrono>
+
+
+#include <bits/stdc++.h>
+
+
+#include "micro.hpp"
+
+
+using namespace std;
+using namespace std::chrono;
+
+
+double eps_vs_t(double time, double t_final) {
+	const double eps_max = 1.0e-1;
+	double eps = eps_max * time;
+	return eps;
+}
+
+
+int main(int argc, char **argv)
+{
+
+	if (argc < 2) {
+		cerr    << "Usage: " << argv[0]
+			<< " n [print=0(def)|1] [steps] [fe_lin=0(def)|fe_one_way|fe_full|rule_mix]" << endl;
+		return 1;
+	}
+
+	const int n = atoi(argv[1]);
+	const int print = (argc > 2 ? atoi(argv[2]) : 0);  // Optional value
+	const int time_steps = (argc > 3 ? atoi(argv[3]) : 10);  // Optional value
+	const int coupling = (argc > 4 ? atoi(argv[4]) : 0);  // Optional value
+
+	const int dir = 0;
+	const double t_final = 0.15;
+	const double dt = t_final / time_steps;
+	double time = 0.0;
+
+	assert(n > 1 && 0 <= print && print < 2 && time_steps >= 0);
+
+	micropp_params_t mic_params;
+
+	mic_params.ngp = 1;
+	mic_params.coupling = new int[1];
+	mic_params.coupling[0] = coupling;
+
+	mic_params.size[0] = n;
+	mic_params.size[1] = n;
+	mic_params.size[2] = n;
+	mic_params.type = MIC3D_8;
+	mic_params.geo_params[0] = 0.1;
+	mic_params.geo_params[1] = 0.02;
+	mic_params.geo_params[2] = 0.01;
+	material_set(&mic_params.materials[0], 0, 1.0e7, 0.3, 0.0, 0.0, 0.0);
+	material_set(&mic_params.materials[1], 0, 3.0e7, 0.3, 0.0, 0.0, 0.0);
+	material_set(&mic_params.materials[2], 0, 3.0e7, 0.3, 0.0, 0.0, 0.0);
+	mic_params.lin_stress = true;
+
+	//mic_params.print();
+
+	micropp<3> micro(mic_params);
+	micro.print_info();
+
+	delete [] mic_params.coupling;
+
+	if (print) {
+		char filename[128];
+		snprintf(filename, 128, "micropp_%d", 0);
+		micro.output (0, filename);
+	}
+
+	ofstream file;
+	file.open("result.dat");
+
+	auto start = high_resolution_clock::now();
+
+	double sig[6];
+	double eps[6] = { 0. };
+
+	cout << scientific;
+
+	for (int t = 0; t < time_steps; ++t) {
+
+		cout << "Time step = " << t << endl;
+
+		eps[dir] = eps_vs_t(time, t_final);
+
+		micro.set_strain(0, eps);
+
+		cout << "eps = ";
+		for (int i = 0; i < 6; ++i) {
+			cout << eps[i] << " ";
+		}
+		cout << endl;
+
+		cout << "Homogenizing ..." << endl;
+		micro.homogenize();
+
+		int non_linear = micro.is_non_linear(0);
+		int cost = micro.get_cost(0);
+		bool has_converged = micro.has_converged(0);
+
+		cout 	<< "NL        = " << non_linear << endl
+			<< "Cost      = " << cost << endl
+			<< "Converged = " << has_converged << endl;
+
+		cout << "sig = ";
+		micro.get_stress(0, sig);
+		for (int i = 0; i < 6; ++i) {
+			cout << sig[i] << "\t";
+		}
+		cout << endl;
+
+		micro.update_vars();
+
+		file    << setw(14)
+			<< eps[dir] << "\t"
+			<< sig[dir] << "\t" << endl;
+
+		if (print) {
+			char filename[128];
+			snprintf(filename, 128, "micropp_%d", t);
+			micro.output (0, filename);
+		}
+
+		time += dt;
+	}
+
+	auto stop = high_resolution_clock::now();
+	auto duration = duration_cast<milliseconds>(stop - start);
+	cout << "time = " << duration.count() << " ms" << endl;
+
+	return 0;
+}