Preparing code to add Rule of Mixtures

include/micro.hpp

@@ -156,10 +156,10 @@ enum {
  */
 
 enum {
-       	LINEAR,
-       	ONE_WAY,
-       	FULL,
-       	RULE_MIXTURE_1
+       	FE_LINEAR,
+       	FE_ONE_WAY,
+       	FE_FULL,
+       	RULE_MIXTURE_LIN_1
 };
 
 
@@ -193,7 +193,7 @@ class micropp {
 		double geo_params[num_geo_params];
 
 		material_t *material_list[MAX_MATERIALS];
-		double ctan_lin[nvoi * nvoi];
+		double ctan_lin_fe[nvoi * nvoi];
 
 		int *elem_type;
 		double *elem_stress;
@@ -217,9 +217,10 @@ class micropp {
 		const bool lin_stress;
 
 		/* Number of micro-problems depending on the type */
-		int num_no_coupling = 0;
-		int num_one_way = 0;
-		int num_full = 0;
+		int num_fe_linear = 0;
+		int num_fe_one_way = 0;
+		int num_fe_full = 0;
+		int num_fe_points = 0;
 
 		/* Linear jacobian for optimization */
 		bool use_A0;
@@ -241,11 +242,18 @@ class micropp {
 
 		/* Private function members */
 
+		/* Linear homogenizations */
 		void homogenize_linear(gp_t<tdim> *gp_ptr);
-		void homogenize_non_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();
+		void calc_ctan_lin_fe_models();
+		void calc_ctan_lin_rule_mixture_lin_1(double ctan[nvoi * nvoi]);
 
 		material_t *get_material(const int e) const;
 

src/homogenize.cpp

@@ -21,9 +21,11 @@
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  */
 
+
 #include <cmath>
 #include <cassert>
 
+
 #include "instrument.hpp"
 #include "micro.hpp"
 
@@ -65,11 +67,7 @@ void micropp<tdim>::homogenize()
 
 		gp_t<tdim> *gp_ptr = &gp_list[igp];
 
-		if (gp_ptr->coupling == RULE_MIXTURE_1) {
-
-			homogenize_rule_mixture_1(gp_ptr);
-
-		} else if (gp_ptr->coupling == LINEAR) {
+		if (gp_ptr->coupling == FE_LINEAR || gp_ptr->coupling == RULE_MIXTURE_LIN_1) {
 
 			/*
 			 * Computational cheap calculation
@@ -78,9 +76,13 @@ void micropp<tdim>::homogenize()
 
 			homogenize_linear(gp_ptr);
 
-		} else {
+		} else if (gp_ptr->coupling == FE_ONE_WAY) {
+
+			homogenize_fe_one_way(gp_ptr);
+
+		} else if (gp_ptr->coupling == FE_FULL) {
 
-			homogenize_non_linear(gp_ptr);
+			homogenize_fe_full(gp_ptr);
 
 		}
 	}
@@ -97,15 +99,98 @@ 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] += gp_ptr->ctan[i * nvoi + j] * gp_ptr->strain[j];
 		}
 	}
 }
 
 
 template<int tdim>
-void micropp<tdim>::homogenize_non_linear(gp_t<tdim> * gp_ptr)
+void micropp<tdim>::homogenize_fe_one_way(gp_t<tdim> * gp_ptr)
+{
+
+	ell_matrix A;  // Jacobian
+	const int ns[3] = { nx, ny, nz };
+	ell_init(&A, dim, dim, ns, CG_ABS_TOL, CG_REL_TOL, CG_MAX_ITS);
+	double *b = (double *) calloc(nndim, sizeof(double));
+	double *du = (double *) calloc(nndim, sizeof(double));
+	double *u = (double *) calloc(nndim, sizeof(double));
+	double *vars_new_aux = (double *) calloc(nvars, sizeof(double));
+
+	double *vars_new = (gp_ptr->allocated) ? gp_ptr->vars_k : vars_new_aux;
+
+	gp_ptr->cost = 0;
+	gp_ptr->subiterated = false;
+
+	// SIGMA 1 Newton-Raphson
+	memcpy(u, gp_ptr->u_n, nndim * sizeof(double));
+
+	newton_t newton = newton_raphson(&A, b, u, du, gp_ptr->strain, gp_ptr->vars_n);
+
+	memcpy(gp_ptr->u_k, u, nndim * sizeof(double));
+	gp_ptr->cost += newton.solver_its;
+	gp_ptr->converged = newton.converged;
+
+	/*
+	 * In case it has not converged do the sub-iterations
+	 */
+	if (gp_ptr->converged == false && subiterations == true) {
+
+		double eps_sub[nvoi], deps_sub[nvoi];
+		memcpy(u, gp_ptr->u_n, nndim * sizeof(double));
+		memcpy(eps_sub, gp_ptr->strain_old, nvoi * sizeof(double));
+		gp_ptr->subiterated = true;
+
+		for (int i = 0; i < nvoi; ++i)
+			deps_sub[i] = (gp_ptr->strain[i] - gp_ptr->strain_old[i]) / nsubiterations;
+
+		for (int its = 0; its < nsubiterations; ++its) {
+
+			for (int j = 0; j < nvoi; ++j) {
+				eps_sub[j] += deps_sub[j];
+			}
+
+			newton = newton_raphson(&A, b, u, du, eps_sub, gp_ptr->vars_n);
+			gp_ptr->cost += newton.solver_its;
+		}
+
+		gp_ptr->converged = newton.converged;
+		memcpy(gp_ptr->u_k, u, nndim * sizeof(double));
+	}
+
+	if (lin_stress) {
+
+		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] += gp_ptr->ctan[i * nvoi + j] * gp_ptr->strain[j];
+			}
+		}
+
+	} else {
+		calc_ave_stress(gp_ptr->u_k, gp_ptr->stress, gp_ptr->vars_n);
+	}
+
+	// Updates <vars_new>
+	bool non_linear = calc_vars_new(gp_ptr->u_k, gp_ptr->vars_n, vars_new);
+
+	if (non_linear == true) {
+		if (gp_ptr->allocated == false) {
+			gp_ptr->allocate();
+			memcpy(gp_ptr->vars_k, vars_new, nvars * sizeof(double));
+		}
+	}
+
+	ell_free(&A);
+	free(b);
+	free(u);
+	free(du);
+	free(vars_new_aux);
+}
+
+
+template<int tdim>
+void micropp<tdim>::homogenize_fe_full(gp_t<tdim> * gp_ptr)
 {
 
 	ell_matrix A;  // Jacobian
@@ -124,8 +209,7 @@ void micropp<tdim>::homogenize_non_linear(gp_t<tdim> * gp_ptr)
 	// SIGMA 1 Newton-Raphson
 	memcpy(u, gp_ptr->u_n, nndim * sizeof(double));
 
-	newton_t newton = newton_raphson(&A, b, u, du, gp_ptr->strain,
-					 gp_ptr->vars_n);
+	newton_t newton = newton_raphson(&A, b, u, du, gp_ptr->strain, gp_ptr->vars_n);
 
 	memcpy(gp_ptr->u_k, u, nndim * sizeof(double));
 	gp_ptr->cost += newton.solver_its;
@@ -162,8 +246,7 @@ void micropp<tdim>::homogenize_non_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] += gp_ptr->ctan[i * nvoi + j] * gp_ptr->strain[j];
 			}
 		}
 
@@ -181,7 +264,7 @@ void micropp<tdim>::homogenize_non_linear(gp_t<tdim> * gp_ptr)
 		}
 	}
 
-	if (gp_ptr->allocated && gp_ptr->coupling == FULL) {
+	if (gp_ptr->allocated) {
 
 		// CTAN 3/6 Newton-Raphsons in 2D/3D
 		double eps_1[6], sig_0[6], sig_1[6];
@@ -201,8 +284,7 @@ void micropp<tdim>::homogenize_non_linear(gp_t<tdim> * gp_ptr)
 			calc_ave_stress(u, sig_1, gp_ptr->vars_n);
 
 			for (int v = 0; v < nvoi; ++v)
-				gp_ptr->ctan[v * nvoi + i] =
-					(sig_1[v] - sig_0[v]) / D_EPS_CTAN_AVE;
+				gp_ptr->ctan[v * nvoi + i] = (sig_1[v] - sig_0[v]) / D_EPS_CTAN_AVE;
 
 		}
 	}

src/micro_common.cpp

@@ -82,27 +82,29 @@ micropp<tdim>::micropp(const micropp_params_t &params):
 		if(params.coupling != nullptr) {
 			gp_list[gp].coupling = params.coupling[gp];
 			switch (params.coupling[gp]) {
-				case LINEAR:
-					num_no_coupling ++;
+				case FE_LINEAR:
+					num_fe_linear ++;
 					break;
-				case ONE_WAY:
-					num_one_way ++;
+				case FE_ONE_WAY:
+					num_fe_one_way ++;
 					break;
-				case FULL:
-					num_full ++;
+				case FE_FULL:
+					num_fe_full ++;
 					break;
 			}
 		} else {
-			gp_list[gp].coupling = ONE_WAY;
-			num_one_way ++;
+			gp_list[gp].coupling = FE_ONE_WAY;
+			num_fe_one_way ++;
 		}
 
+		num_fe_points = num_fe_linear + num_fe_one_way + num_fe_full;
+
 		gp_list[gp].nndim = nndim;
 		gp_list[gp].nvars = nvars;
 
 		if (params.coupling == nullptr ||
-		    (params.coupling[gp] == ONE_WAY ||
-		     params.coupling[gp] == FULL)) {
+		    (params.coupling[gp] == FE_ONE_WAY ||
+		     params.coupling[gp] == FE_FULL)) {
 			gp_list[gp].allocate_u();
 		}
 	}
@@ -140,8 +142,7 @@ micropp<tdim>::micropp(const micropp_params_t &params):
 
 #pragma omp parallel for schedule(dynamic,1)
 		for (int i = 0; i < num_of_A0s; ++i) {
-			ell_init(&A0[i], dim, dim, params.size, CG_ABS_TOL,
-				 CG_REL_TOL, CG_MAX_ITS);
+			ell_init(&A0[i], dim, dim, params.size, CG_ABS_TOL, CG_REL_TOL, CG_MAX_ITS);
 			double *u = (double *) calloc(nndim, sizeof(double));
 			assembly_mat(&A0[i], u, nullptr);
 			free(u);
@@ -150,14 +151,27 @@ micropp<tdim>::micropp(const micropp_params_t &params):
 
 	/* Average tangent constitutive tensor initialization */
 
-	memset(ctan_lin, 0.0, nvoi * nvoi * sizeof(double));
+	memset(ctan_lin_fe, 0.0, nvoi * nvoi * sizeof(double));
 
 	if (calc_ctan_lin_flag) {
-		calc_ctan_lin();
+		if (num_fe_points > 0) {
+			calc_ctan_lin_fe_models();
+		}
 	}
 
 	for (int gp = 0; gp < ngp; ++gp) {
-		memcpy(gp_list[gp].ctan, ctan_lin, nvoi * nvoi * sizeof(double));
+		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) {
+
+			double ctan[nvoi * nvoi];
+			calc_ctan_lin_rule_mixture_lin_1(ctan);
+			memcpy(gp_list[gp].ctan, ctan, nvoi * nvoi * sizeof(double));
+
+		}
 	}
 
 	/* Open the log file */
@@ -170,9 +184,7 @@ micropp<tdim>::micropp(const micropp_params_t &params):
 		strcpy(filename, file_name_string.c_str());
 
 		ofstream_log.open(filename, ios::out);
-		ofstream_log
-			<< "#<gp_id>  <non-linear>  <cost>  <converged>"
-			<< endl;
+		ofstream_log << "#<gp_id>  <non-linear>  <cost>  <converged>" << endl;
 	}
 
 }
@@ -257,7 +269,7 @@ int micropp<tdim>::get_non_linear_gps(void) const
 
 
 template <int tdim>
-void micropp<tdim>::calc_ctan_lin()
+void micropp<tdim>::calc_ctan_lin_fe_models()
 {
 
 #pragma omp parallel for schedule(dynamic,1)
@@ -280,7 +292,7 @@ void micropp<tdim>::calc_ctan_lin()
 		calc_ave_stress(u, sig);
 
 		for (int v = 0; v < nvoi; ++v) {
-			ctan_lin[v * nvoi + i] = sig[v] / D_EPS_CTAN_AVE;
+			ctan_lin_fe[v * nvoi + i] = sig[v] / D_EPS_CTAN_AVE;
 		}
 
 		ell_free(&A);
@@ -291,6 +303,12 @@ void micropp<tdim>::calc_ctan_lin()
 }
 
 
+template <int tdim>
+void micropp<tdim>::calc_ctan_lin_rule_mixture_lin_1(double ctan[nvoi * nvoi])
+{
+}
+
+
 template <int tdim>
 material_t *micropp<tdim>::get_material(const int e) const
 {
@@ -819,15 +837,17 @@ 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 << "ONE_WAY     : " << num_one_way     << " GPs" << endl;
-	cout << "FULL        : " << num_full        << " GPs" << endl;
-	cout << "USE A0      : " << use_A0 << endl;
-	cout << "NUM SUBITS  : " << nsubiterations << endl;
-	cout << "MPI RANK    : " << mpi_rank << endl;
+
+	cout << "MATRIX [%]        : " << Vm << endl;
+	cout << "FIBER  [%]        : " << Vf << endl;
+
+	cout << "FE_LINEAR         : " << num_fe_linear    << " GPs" << endl;
+	cout << "FE_ONE_WAY        : " << num_fe_one_way   << " GPs" << endl;
+	cout << "FE_FULL           : " << num_fe_full      << " GPs" << endl;
+	cout << "USE A0            : " << use_A0 << endl;
+	cout << "NUM SUBITS        : " << nsubiterations << endl;
+	cout << "MPI RANK          : " << mpi_rank << endl;
+
 #ifdef _OPENACC
 	int acc_num_gpus = acc_get_num_devices(acc_device_nvidia);
 	cout << "ACC NUM GPUS      : " << acc_num_gpus << endl;
@@ -856,10 +876,10 @@ void micropp<tdim>::print_info() const
 	}
 
 	cout << endl;
-	cout << "ctan lin = " << endl;
+	cout << "ctan_lin_fe = " << endl;
 	for (int i = 0; i < 6; ++i) {
 		for (int j = 0; j < 6; ++j) {
-			cout << ctan_lin[i * 6 + j] << "\t";
+			cout << ctan_lin_fe[i * 6 + j] << "\t";
 		}
 		cout << endl;
 	}