16 #include "../general/forall.hpp" 82 Vector products_no_m(m_cnt); products_no_m = 1.0;
83 for (
int m_p = 0; m_p < m_cnt; m_p++)
85 for (
int m_a = 0; m_a < m_cnt; m_a++)
87 if (m_p != m_a) { products_no_m(m_p) *= averages(m_a); }
90 const double pnm_sum = products_no_m.
Sum();
92 if (pnm_sum == 0.0) { weights = 1.0 / m_cnt;
return; }
93 for (
int m = 0; m < m_cnt; m++) { weights(m) = products_no_m(m) / pnm_sum; }
95 MFEM_ASSERT(fabs(weights.
Sum() - 1.0) < 1e-14,
96 "Error: sum should be 1 always: " << weights.
Sum());
113 for (
int e = 0; e < NE; e++)
132 for (
int q = 0; q < nsp; q++)
142 const double w_detA = ip.
weight * A.
Det();
143 for (
int m = 0; m < m_cnt; m++)
146 averages(m) +=
tmop_q_arr[m]->EvalW(T) * w_detA;
157 MPI_Allreduce(MPI_IN_PLACE, averages.
GetData(), m_cnt,
158 MPI_DOUBLE, MPI_SUM, par_nodes->ParFESpace()->GetComm());
159 MPI_Allreduce(MPI_IN_PLACE, &volume, 1, MPI_DOUBLE, MPI_SUM,
160 par_nodes->ParFESpace()->GetComm());
171 double metric = metric_tilde;
174 metric = std::pow(metric_tilde,
exponent);
179 metric = metric_tilde/(beta-metric_tilde);
187 double denominator = 1.0;
194 double detT = Jpt.
Det();
224 MFEM_VERIFY(
Jtr != NULL,
225 "Requires a target Jacobian, use SetTargetJacobian().");
234 const double cos_Jpr = (col1 * col2) / norm_prod,
235 sin_Jpr = fabs(Jpr.
Det()) / norm_prod;
240 const double cos_Jtr = (col1 * col2) / norm_prod,
241 sin_Jtr = fabs(
Jtr->
Det()) / norm_prod;
243 return 0.5 * (1.0 - cos_Jpr * cos_Jtr - sin_Jpr * sin_Jtr);
248 MFEM_VERIFY(
Jtr != NULL,
249 "Requires a target Jacobian, use SetTargetJacobian().");
258 double norm_c1 = col1.
Norml2(),
261 double cos_Jpr_12 = (col1 * col2) / (norm_c1 * norm_c2),
262 cos_Jpr_13 = (col1 * col3) / (norm_c1 * norm_c3),
263 cos_Jpr_23 = (col2 * col3) / (norm_c2 * norm_c3);
264 double sin_Jpr_12 = std::sqrt(1.0 - cos_Jpr_12 * cos_Jpr_12),
265 sin_Jpr_13 = std::sqrt(1.0 - cos_Jpr_13 * cos_Jpr_13),
266 sin_Jpr_23 = std::sqrt(1.0 - cos_Jpr_23 * cos_Jpr_23);
274 double cos_Jtr_12 = (col1 * col2) / (norm_c1 * norm_c2),
275 cos_Jtr_13 = (col1 * col3) / (norm_c1 * norm_c3),
276 cos_Jtr_23 = (col2 * col3) / (norm_c2 * norm_c3);
277 double sin_Jtr_12 = std::sqrt(1.0 - cos_Jtr_12 * cos_Jtr_12),
278 sin_Jtr_13 = std::sqrt(1.0 - cos_Jtr_13 * cos_Jtr_13),
279 sin_Jtr_23 = std::sqrt(1.0 - cos_Jtr_23 * cos_Jtr_23);
281 return (3.0 - cos_Jpr_12 * cos_Jtr_12 - sin_Jpr_12 * sin_Jtr_12
282 - cos_Jpr_13 * cos_Jtr_13 - sin_Jpr_13 * sin_Jtr_13
283 - cos_Jpr_23 * cos_Jtr_23 - sin_Jpr_23 * sin_Jtr_23) / 6.0;
288 MFEM_VERIFY(
Jtr != NULL,
289 "Requires a target Jacobian, use SetTargetJacobian().");
303 return 0.5 * (ratio_Jpr / ratio_Jtr + ratio_Jtr / ratio_Jpr) - 1.0;
308 MFEM_VERIFY(
Jtr != NULL,
309 "Requires a target Jacobian, use SetTargetJacobian().");
318 double norm_c1 = col1.
Norml2(),
321 double ratio_Jpr_1 = norm_c1 / std::sqrt(norm_c2 * norm_c3),
322 ratio_Jpr_2 = norm_c2 / std::sqrt(norm_c1 * norm_c3),
323 ratio_Jpr_3 = norm_c3 / std::sqrt(norm_c1 * norm_c2);
331 double ratio_Jtr_1 = norm_c1 / std::sqrt(norm_c2 * norm_c3),
332 ratio_Jtr_2 = norm_c2 / std::sqrt(norm_c1 * norm_c3),
333 ratio_Jtr_3 = norm_c3 / std::sqrt(norm_c1 * norm_c2);
335 return ( 0.5 * (ratio_Jpr_1 / ratio_Jtr_1 + ratio_Jtr_1 / ratio_Jpr_1) +
336 0.5 * (ratio_Jpr_2 / ratio_Jtr_2 + ratio_Jtr_2 / ratio_Jpr_2) +
337 0.5 * (ratio_Jpr_3 / ratio_Jtr_3 + ratio_Jtr_3 / ratio_Jpr_3) - 3.0
343 return 0.5 * Jpt.
FNorm2() / Jpt.
Det() - 1.0;
424 const double c2 = weight*c1*c1;
468 MFEM_VERIFY(
Jtr != NULL,
469 "Requires a target Jacobian, use SetTargetJacobian().");
473 Id(0,0) = 1;
Id(0,1) = 0;
474 Id(1,0) = 0;
Id(1,1) = 1;
530 const double c2 = weight*c1/2;
531 const double c3 = c1*c2;
547 return 0.5*I1b*I1b - 2.0;
605 return 0.5*(I2b + 1.0/I2b) - 1.0;
636 double det = Jpt.
Det();
638 return JtJ.
FNorm2()/(det*det) - 2*Jpt.
FNorm2()/det + 2.0;
646 return I1b*(I1b - 2.0);
674 return 0.5*(I2b*I2b + 1./(I2b*I2b) - 2.);
693 const double I2 =
ie.
Get_I2(), I2inv_sq = 1.0 / (I2 * I2);
701 MFEM_VERIFY(
Jtr != NULL,
702 "Requires a target Jacobian, use SetTargetJacobian().");
708 Id(0,0) = 1;
Id(0,1) = 0;
709 Id(1,0) = 0;
Id(1,1) = 1;
710 Id *=
Mat.FNorm()/pow(2,0.5);
719 MFEM_VERIFY(
Jtr != NULL,
720 "Requires a target Jacobian, use SetTargetJacobian().");
724 Id(0,0) = 1;
Id(0,1) = 0;
725 Id(1,0) = 0;
Id(1,1) = 1;
739 return (I2b - 1.0)*(I2b - 1.0) - I2b + std::sqrt(I2b*I2b +
eps);
744 MFEM_ABORT(
"Metric not implemented yet. Use metric mu_55 instead.");
752 MFEM_ABORT(
"Metric not implemented yet. Use metric mu_55 instead.");
760 return 0.5*(I2b - 1.0)*(I2b - 1.0)/(I2b -
tau0);
771 const double c = (I2b - 1.0)/(I2b -
tau0);
787 const double c0 = 1.0/(I2b -
tau0);
788 const double c = c0*(I2b - 1.0);
799 return Jpt.FNorm() * inv.FNorm() / 3.0 - 1.0;
847 const double a = weight/(6*std::sqrt(I1b_I2b));
890 const double c1 = weight/9;
900 return Jpt.
FNorm2() / 3.0 / pow(Jpt.
Det(), 2.0 / 3.0) - 1.0;
933 const double fnorm = Jpt.FNorm();
934 return fnorm * fnorm * fnorm / pow(3.0, 1.5) / Jpt.
Det() - 1.0;
970 return (I3b - 1.0)*(I3b - 1.0) - I3b + std::sqrt(I3b*I3b +
eps);
977 const double c = 2*I3b-3+(I3b)/(std::pow((I3b*I3b+
eps),0.5));
989 const double c0 = I3b*I3b+
eps;
990 const double c1 = 2 + 1/(pow(c0,0.5)) - I3b*I3b/(pow(c0,1.5));
991 const double c2 = 2*I3b - 3 + I3b/(pow(c0,0.5));
1011 const double c = std::pow(d, -2.0/3.0);
1018 MFEM_ABORT(
"Metric not implemented yet.");
1023 const double weight,
1026 MFEM_ABORT(
"Metric not implemented yet.");
1047 const double weight,
1061 return 0.5 * (Jpt.
Det() + 1.0 / Jpt.
Det()) - 1.0;
1069 return 0.5*(I3b + 1.0/I3b) - 1.0;
1082 const double weight,
1100 invt.
Add(-1.0, Jpt);
1128 const double weight,
1142 const double c1 = weight*c0*c0;
1143 const double c2 = -2*c0*c1;
1158 adj_J_t.
Add(1.0, Jpt);
1159 return 1.0 / 6.0 / Jpt.
Det() * adj_J_t.
FNorm2();
1204 const double p13 = weight * pow(
ie.
Get_I3b(), 1.0/3.0),
1205 m13 = weight * pow(
ie.
Get_I3b(), -1.0/3.0),
1206 m23 = weight * pow(
ie.
Get_I3b(), -2.0/3.0),
1207 m43 = weight * pow(
ie.
Get_I3b(), -4.0/3.0),
1208 m53 = weight * pow(
ie.
Get_I3b(), -5.0/3.0),
1209 m73 = weight * pow(
ie.
Get_I3b(), -7.0/3.0);
1229 double fnorm = Jpt.FNorm();
1230 return fnorm * fnorm * fnorm - 3.0 * sqrt(3.0) * (log(Jpt.
Det()) + 1.0);
1270 return 0.5*(I3b - 1.0)*(I3b - 1.0)/(I3b -
tau0);
1281 const double c = (I3b - 1.0)/(I3b -
tau0);
1287 const double weight,
1302 const double c0 = 1.0/(I3b -
tau0);
1303 const double c = c0*(I3b - 1.0);
1314 const double fnorm = Jpt.FNorm();
1315 return fnorm * fnorm * fnorm / pow(3.0, 1.5) - Jpt.
Det();
1349 MFEM_VERIFY(
Jtr != NULL,
1350 "Requires a target Jacobian, use SetTargetJacobian().");
1365 Mult(AdjAt, WtW, WRK);
1375 MFEM_VERIFY(
Jtr != NULL,
1376 "Requires a target Jacobian, use SetTargetJacobian().");
1383 double sqalpha = pow(Jpr.
Det(), 0.5),
1384 sqomega = pow(
Jtr->
Det(), 0.5);
1386 return 0.5*pow(sqalpha/sqomega - sqomega/sqalpha, 2.);
1391 MFEM_VERIFY(
Jtr != NULL,
1392 "Requires a target Jacobian, use SetTargetJacobian().");
1407 MFEM_VERIFY(
Jtr != NULL,
1408 "Requires a target Jacobian, use SetTargetJacobian().");
1416 aw = Jpr.FNorm()/
Jtr->FNorm();
1428 MFEM_VERIFY(
nodes,
"Nodes are not given!");
1429 MFEM_ASSERT(
avg_volume == 0.0,
"The average volume is already computed!");
1432 const int NE = mesh->
GetNE();
1434 double volume = 0.0;
1436 for (
int i = 0; i < NE; i++)
1460 area_NE[0] = volume; area_NE[1] = NE;
1461 MPI_Allreduce(area_NE, area_NE + 2, 2, MPI_DOUBLE, MPI_SUM,
comm);
1480 const int NE = mesh->
GetNE();
1482 if (NE == 0) {
return; }
1485 "mixed meshes are not supported");
1486 MFEM_VERIFY(!fes.
IsVariableOrder(),
"variable orders are not supported");
1488 const int sdim = fes.
GetVDim();
1489 const int nvdofs = sdim*fe.GetDof();
1491 xe.
Size() == NE*nvdofs,
"invalid input Vector 'xe'!");
1501 if (dof_map->
Size() == 0) { dof_map =
nullptr; }
1505 Vector elfun_lex, elfun_nat;
1513 for (
int e = 0; e < NE; e++)
1524 const int ndofs = fe.GetDof();
1525 for (
int d = 0; d < sdim; d++)
1527 for (
int i_lex = 0; i_lex < ndofs; i_lex++)
1529 elfun_nat[(*dof_map)[i_lex]+d*ndofs] =
1530 elfun_lex[i_lex+d*ndofs];
1549 default: MFEM_ABORT(
"TargetType not added to ContainsVolumeInfo.");
1559 MFEM_CONTRACT_VAR(elfun);
1567 MFEM_ASSERT(Wideal.
Width() == Jtr.
SizeJ(),
"");
1573 for (
int i = 0; i < ir.
GetNPoints(); i++) { Jtr(i) = Wideal; }
1585 el_volume =
avg_volume / ncmesh->GetElementSizeReduction(e_id);
1589 1./W.Height()), Wideal);
1590 for (
int i = 0; i < ir.
GetNPoints(); i++) { Jtr(i) = W; }
1613 const double det = Jtr(i).Det();
1614 MFEM_VERIFY(det > 0.0,
"The given mesh is inverted!");
1615 Jtr(i).Set(std::pow(det / detW, 1./
dim), Wideal);
1621 MFEM_ABORT(
"invalid target type!");
1631 MFEM_CONTRACT_VAR(elfun);
1661 "Target type GIVEN_FULL requires a MatrixCoefficient.");
1678 MFEM_ABORT(
"Incompatible target type for analytic adaptation!");
1696 "Target type GIVEN_FULL requires a TMOPMatrixCoefficient.");
1698 for (
int d = 0; d < fe->
GetDim(); d++)
1711 MFEM_ABORT(
"Incompatible target type for analytic adaptation!");
1721 static inline void device_copy(
double *d_dest,
const double *d_src,
int size)
1723 MFEM_FORALL(i, size, d_dest[i] = d_src[i];);
1731 MFEM_VERIFY(
adapt_eval,
"SetAdaptivityEvaluator() has not been called!")
1732 MFEM_VERIFY(
ncomp > 0,
"No target specifications have been set!");
1776 MFEM_VERIFY(ndof ==
tspec.
Size()/
ncomp,
"Inconsistency in SetTspecAtIndex.");
1778 const auto tspec__d = tspec_.
Read();
1780 const int offset = idx*ndof;
1781 internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1788 "Discrete target size should be ordered byNodes.");
1798 "Discrete target skewness should be ordered byNodes.");
1808 "Discrete target aspect ratio should be ordered byNodes.");
1818 "Discrete target orientation should be ordered byNodes.");
1829 #endif // MFEM_USE_MPI 1846 const auto tspec_temp_d = tspec_temp.
Read();
1848 internal::device_copy(tspec_d, tspec_temp_d, tspec_temp.
Size());
1850 const auto tspec__d = tspec_.
Read();
1851 const int offset = (
ncomp-vdim)*ndof;
1852 internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1859 MFEM_VERIFY(ndof ==
tspec.
Size()/
ncomp,
"Inconsistency in SetTargetSpec.");
1861 const auto tspec__d = tspec_.
Read();
1863 const int offset = idx*ndof;
1864 internal::device_copy(tspec_d + offset, tspec__d, ndof*vdim);
1871 "Discrete target size should be ordered byNodes.");
1881 "Discrete target skewness should be ordered byNodes.");
1891 "Discrete target aspect ratio should be ordered byNodes.");
1901 "Discrete target orientation should be ordered byNodes.");
1910 MFEM_VERIFY(
adapt_eval,
"SetAdaptivityEvaluator() has not been called!")
1911 MFEM_VERIFY(
ncomp > 0,
"No target specifications have been set!");
1931 if (idx < 0) {
return; }
1935 "Inconsistency in GetSerialDiscreteTargetSpec.");
1937 for (
int i = 0; i < ndof*vdim; i++)
1939 tspec_(i) =
tspec(i + idx*ndof);
1966 MFEM_VERIFY(
tspec.
Size() > 0,
"Target specification is not set!");
1982 int dofidx,
int dir,
1985 MFEM_VERIFY(
tspec.
Size() > 0,
"Target specification is not set!");
1991 for (
int i = 0; i <
ncomp; i++)
1993 tspec(dofs[dofidx]+i*cnt) = IntData(dofs[dofidx] + i*cnt + dir*cnt*
ncomp);
2000 MFEM_VERIFY(
tspec.
Size() > 0,
"Target specification is not set!");
2005 for (
int i = 0; i <
ncomp; i++)
2020 ntspec_dofs = ndofs*
ncomp;
2022 Vector tspec_vals(ntspec_dofs);
2033 MFEM_ABORT(
"Incompatible target type for discrete adaptation!");
2050 MFEM_VERIFY(
tspec_fesv,
"No target specifications have been set.");
2065 ntspec_dofs = ndofs*
ncomp;
2067 Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
2068 par_vals_c1, par_vals_c2, par_vals_c3;
2077 MFEM_VERIFY(
amr_el >= 0,
" Target being constructed for an AMR element.");
2078 for (
int i = 0; i <
ncomp; i++)
2080 for (
int j = 0; j < ndofs; j++)
2094 for (
int q = 0; q < nqp; q++)
2099 for (
int d = 0; d < 4; d++)
2108 double min_size = par_vals.
Min();
2115 MFEM_VERIFY(min_size > 0.0,
2116 "Non-positive size propagated in the target definition.");
2118 const double size = std::max(shape * par_vals, min_size);
2119 Jtr(q).Set(std::pow(size, 1.0/
dim), Jtr(q));
2132 const double min_size = par_vals.
Min();
2133 MFEM_VERIFY(min_size > 0.0,
2134 "Non-positive aspect-ratio propagated in the target definition.");
2136 const double aspectratio = shape * par_vals;
2138 D_rho(0,0) = 1./pow(aspectratio,0.5);
2139 D_rho(1,1) = pow(aspectratio,0.5);
2149 const double rho1 = shape * par_vals_c1;
2150 const double rho2 = shape * par_vals_c2;
2151 const double rho3 = shape * par_vals_c3;
2153 D_rho(0,0) = pow(rho1,2./3.);
2154 D_rho(1,1) = pow(rho2,2./3.);
2155 D_rho(2,2) = pow(rho3,2./3.);
2160 Mult(D_rho, Temp, Jtr(q));
2170 const double skew = shape * par_vals;
2174 Q_phi(0,1) = cos(skew);
2175 Q_phi(1,1) = sin(skew);
2185 const double phi12 = shape * par_vals_c1;
2186 const double phi13 = shape * par_vals_c2;
2187 const double chi = shape * par_vals_c3;
2191 Q_phi(0,1) = cos(phi12);
2192 Q_phi(0,2) = cos(phi13);
2194 Q_phi(1,1) = sin(phi12);
2195 Q_phi(1,2) = sin(phi13)*cos(chi);
2197 Q_phi(2,2) = sin(phi13)*sin(chi);
2202 Mult(Q_phi, Temp, Jtr(q));
2212 const double theta = shape * par_vals;
2213 R_theta(0,0) = cos(theta);
2214 R_theta(0,1) = -sin(theta);
2215 R_theta(1,0) = sin(theta);
2216 R_theta(1,1) = cos(theta);
2226 const double theta = shape * par_vals_c1;
2227 const double psi = shape * par_vals_c2;
2228 const double beta = shape * par_vals_c3;
2230 double ct = cos(theta), st = sin(theta),
2231 cp = cos(psi), sp = sin(psi),
2232 cb = cos(beta), sb = sin(beta);
2235 R_theta(0,0) = ct*sp;
2236 R_theta(1,0) = st*sp;
2239 R_theta(0,1) = -st*cb + ct*cp*sb;
2240 R_theta(1,1) = ct*cb + st*cp*sb;
2241 R_theta(2,1) = -sp*sb;
2243 R_theta(0,0) = -st*sb - ct*cp*cb;
2244 R_theta(1,0) = ct*sb - st*cp*cb;
2245 R_theta(2,0) = sp*cb;
2248 Jtrcomp_q = R_theta;
2250 Mult(R_theta, Temp, Jtr(q));
2256 MFEM_ABORT(
"Incompatible target type for discrete adaptation!");
2267 MFEM_VERIFY(
tspec_fesv,
"No target specifications have been set.");
2282 ntspec_dofs = ndofs*
ncomp;
2284 Vector shape(ndofs), tspec_vals(ntspec_dofs), par_vals,
2285 par_vals_c1(ndofs), par_vals_c2(ndofs), par_vals_c3(ndofs);
2295 grad_e_c2(ndofs,
dim),
2296 grad_e_c3(ndofs,
dim);
2297 Vector grad_ptr_c1(grad_e_c1.GetData(), ndofs*
dim),
2298 grad_ptr_c2(grad_e_c2.GetData(), ndofs*
dim),
2317 grad_phys.Mult(par_vals, grad_ptr_c1);
2320 grad_e_c1.MultTranspose(shape, grad_q);
2322 const double min_size = par_vals.
Min();
2323 MFEM_VERIFY(min_size > 0.0,
2324 "Non-positive size propagated in the target definition.");
2325 const double size = std::max(shape * par_vals, min_size);
2326 double dz_dsize = (1./
dim)*pow(size, 1./
dim - 1.);
2328 Mult(Jtrcomp_q, Jtrcomp_d, work1);
2329 Mult(Jtrcomp_r, work1, work2);
2331 for (
int d = 0; d <
dim; d++)
2335 work1.
Set(dz_dsize, work1);
2337 AddMult(work1, work2, dJtr_i);
2350 grad_phys.Mult(par_vals, grad_ptr_c1);
2353 grad_e_c1.MultTranspose(shape, grad_q);
2355 const double aspectratio = shape * par_vals;
2357 dD_rho(0,0) = -0.5*pow(aspectratio,-1.5);
2358 dD_rho(1,1) = 0.5*pow(aspectratio,-0.5);
2360 Mult(Jtrcomp_s, Jtrcomp_r, work1);
2361 Mult(work1, Jtrcomp_q, work2);
2363 for (
int d = 0; d <
dim; d++)
2368 AddMult(work2, work1, dJtr_i);
2375 par_vals_c1.SetData(par_vals.
GetData());
2376 par_vals_c2.SetData(par_vals.
GetData()+ndofs);
2379 grad_phys.Mult(par_vals_c1, grad_ptr_c1);
2380 grad_phys.Mult(par_vals_c2, grad_ptr_c2);
2381 grad_phys.Mult(par_vals_c3, grad_ptr_c3);
2384 grad_e_c1.MultTranspose(shape, grad_q1);
2385 grad_e_c2.MultTranspose(shape, grad_q2);
2388 const double rho1 = shape * par_vals_c1;
2389 const double rho2 = shape * par_vals_c2;
2390 const double rho3 = shape * par_vals_c3;
2392 dD_rho(0,0) = (2./3.)*pow(rho1,-1./3.);
2393 dD_rho(1,1) = (2./3.)*pow(rho2,-1./3.);
2394 dD_rho(2,2) = (2./3.)*pow(rho3,-1./3.);
2396 Mult(Jtrcomp_s, Jtrcomp_r, work1);
2397 Mult(work1, Jtrcomp_q, work2);
2400 for (
int d = 0; d <
dim; d++)
2404 work1(0,0) *= grad_q1(d);
2405 work1(1,2) *= grad_q2(d);
2406 work1(2,2) *= grad_q3(d);
2408 AddMult(work2, work1, dJtr_i);
2420 grad_phys.Mult(par_vals, grad_ptr_c1);
2423 grad_e_c1.MultTranspose(shape, grad_q);
2425 const double skew = shape * par_vals;
2429 dQ_phi(0,1) = -sin(skew);
2430 dQ_phi(1,1) = cos(skew);
2432 Mult(Jtrcomp_s, Jtrcomp_r, work2);
2434 for (
int d = 0; d <
dim; d++)
2439 Mult(work1, Jtrcomp_d, work3);
2440 AddMult(work2, work3, dJtr_i);
2447 par_vals_c1.SetData(par_vals.
GetData());
2448 par_vals_c2.SetData(par_vals.
GetData()+ndofs);
2451 grad_phys.Mult(par_vals_c1, grad_ptr_c1);
2452 grad_phys.Mult(par_vals_c2, grad_ptr_c2);
2453 grad_phys.Mult(par_vals_c3, grad_ptr_c3);
2456 grad_e_c1.MultTranspose(shape, grad_q1);
2457 grad_e_c2.MultTranspose(shape, grad_q2);
2460 const double phi12 = shape * par_vals_c1;
2461 const double phi13 = shape * par_vals_c2;
2462 const double chi = shape * par_vals_c3;
2466 dQ_phi(0,1) = -sin(phi12);
2467 dQ_phi(1,1) = cos(phi12);
2470 dQ_phi13(0,2) = -sin(phi13);
2471 dQ_phi13(1,2) = cos(phi13)*cos(chi);
2472 dQ_phi13(2,2) = cos(phi13)*sin(chi);
2475 dQ_phichi(1,2) = -sin(phi13)*sin(chi);
2476 dQ_phichi(2,2) = sin(phi13)*cos(chi);
2478 Mult(Jtrcomp_s, Jtrcomp_r, work2);
2480 for (
int d = 0; d <
dim; d++)
2484 work1 *= grad_q1(d);
2485 work1.Add(grad_q2(d), dQ_phi13);
2486 work1.Add(grad_q3(d), dQ_phichi);
2487 Mult(work1, Jtrcomp_d, work3);
2488 AddMult(work2, work3, dJtr_i);
2500 grad_phys.Mult(par_vals, grad_ptr_c1);
2503 grad_e_c1.MultTranspose(shape, grad_q);
2505 const double theta = shape * par_vals;
2506 dR_theta(0,0) = -sin(theta);
2507 dR_theta(0,1) = -cos(theta);
2508 dR_theta(1,0) = cos(theta);
2509 dR_theta(1,1) = -sin(theta);
2511 Mult(Jtrcomp_q, Jtrcomp_d, work1);
2512 Mult(Jtrcomp_s, work1, work2);
2513 for (
int d = 0; d <
dim; d++)
2518 AddMult(work1, work2, dJtr_i);
2525 par_vals_c1.SetData(par_vals.
GetData());
2526 par_vals_c2.SetData(par_vals.
GetData()+ndofs);
2529 grad_phys.Mult(par_vals_c1, grad_ptr_c1);
2530 grad_phys.Mult(par_vals_c2, grad_ptr_c2);
2531 grad_phys.Mult(par_vals_c3, grad_ptr_c3);
2534 grad_e_c1.MultTranspose(shape, grad_q1);
2535 grad_e_c2.MultTranspose(shape, grad_q2);
2538 const double theta = shape * par_vals_c1;
2539 const double psi = shape * par_vals_c2;
2540 const double beta = shape * par_vals_c3;
2542 const double ct = cos(theta), st = sin(theta),
2543 cp = cos(psi), sp = sin(psi),
2544 cb = cos(beta), sb = sin(beta);
2547 dR_theta(0,0) = -st*sp;
2548 dR_theta(1,0) = ct*sp;
2551 dR_theta(0,1) = -ct*cb - st*cp*sb;
2552 dR_theta(1,1) = -st*cb + ct*cp*sb;
2555 dR_theta(0,0) = -ct*sb + st*cp*cb;
2556 dR_theta(1,0) = -st*sb - ct*cp*cb;
2564 dR_beta(0,1) = st*sb + ct*cp*cb;
2565 dR_beta(1,1) = -ct*sb + st*cp*cb;
2566 dR_beta(2,1) = -sp*cb;
2568 dR_beta(0,0) = -st*cb + ct*cp*sb;
2569 dR_beta(1,0) = ct*cb + st*cp*sb;
2573 dR_psi(0,0) = ct*cp;
2574 dR_psi(1,0) = st*cp;
2577 dR_psi(0,1) = 0. - ct*sp*sb;
2578 dR_psi(1,1) = 0. + st*sp*sb;
2579 dR_psi(2,1) = -cp*sb;
2581 dR_psi(0,0) = 0. + ct*sp*cb;
2582 dR_psi(1,0) = 0. + st*sp*cb;
2583 dR_psi(2,0) = cp*cb;
2585 Mult(Jtrcomp_q, Jtrcomp_d, work1);
2586 Mult(Jtrcomp_s, work1, work2);
2587 for (
int d = 0; d <
dim; d++)
2591 work1 *= grad_q1(d);
2592 work1.Add(grad_q2(d), dR_psi);
2593 work1.Add(grad_q3(d), dR_beta);
2594 AddMult(work1, work2, dJtr_i);
2602 MFEM_ABORT(
"Incompatible target type for discrete adaptation!");
2621 for (
int j = 0; j <
dim; j++)
2623 for (
int i = 0; i < cnt; i++)
2632 for (
int i = 0; i < cnt; i++)
2643 double dx,
bool use_flag,
2651 totmix = 1+2*(
dim-2);
2660 for (
int j = 0; j <
dim; j++)
2662 for (
int i = 0; i < cnt; i++)
2671 for (
int i = 0; i < cnt; i++)
2680 for (
int k1 = 0; k1 <
dim; k1++)
2682 for (
int k2 = 0; (k1 != k2) && (k2 <
dim); k2++)
2684 for (
int i = 0; i < cnt; i++)
2695 for (
int i = 0; i < cnt; i++)
2726 f.GetVDim(),
f.GetOrdering());
2737 f.GetVDim(),
f.GetOrdering());
2764 PA.H.GetMemory().DeleteDevice(copy_to_host);
2765 PA.H0.GetMemory().DeleteDevice(copy_to_host);
2766 if (!copy_to_host && !
PA.Jtr.GetMemory().HostIsValid())
2768 PA.Jtr_needs_update =
true;
2770 PA.Jtr.GetMemory().DeleteDevice(copy_to_host);
2779 for (
int i = 0; i <
ElemDer.Size(); i++)
2794 "'dist' must be a scalar GridFunction!");
2884 MFEM_VERIFY(
surf_fit_gf,
"Surface fitting has not been enabled.");
2887 double loc_max = 0.0, loc_sum = 0.0;
2893 loc_max = std::max(loc_max, std::abs((*
surf_fit_gf)(i)));
2902 MPI_Allreduce(&loc_max, &err_max, 1, MPI_DOUBLE, MPI_MAX, comm);
2903 MPI_Allreduce(&loc_cnt, &glob_cnt, 1, MPI_INT, MPI_SUM, comm);
2904 MPI_Allreduce(&loc_sum, &err_avg, 1, MPI_DOUBLE, MPI_SUM, comm);
2905 err_avg = err_avg / glob_cnt;
2907 err_avg = loc_sum / loc_cnt;
2943 const int el_id = T.ElementNo;
3007 Vector adapt_lim_gf_q, adapt_lim_gf0_q;
3008 if (adaptive_limiting)
3021 const double weight = ip.
weight * Jtr_i.
Det();
3041 if (adaptive_limiting)
3043 const double diff = adapt_lim_gf_q(i) - adapt_lim_gf0_q(i);
3047 energy += weight * val;
3059 for (
int s = 0;
s < dofs.
Size();
s++)
3066 sigma_e(
s) * sigma_e(
s);
3081 NEsplit = elfun.
Size() / (dof*
dim), el_id = T.ElementNo;
3090 for (
int e = 0; e < NEsplit; e++)
3097 for (
int i = 0; i < dof; i++)
3099 for (
int d = 0; d <
dim; d++)
3105 elfun_child(i + d*dof) = elfun(i + e*dof + d*dof*NEsplit);
3112 double el_energy = 0;
3140 const double weight = ip.
weight * Jtr_i.
Det();
3149 el_energy += weight * val;
3152 energy += el_energy;
3199 const double weight = ip.
weight * Jtr_i.
Det();
3208 energy += weight * val;
3272 Vector shape,
p, p0, d_vals, grad;
3289 d_vals.
SetSize(nqp); d_vals = 1.0;
3313 for (
int q = 0; q < nqp; q++)
3340 for (
int d = 0; d <
dim; d++)
3344 d_detW_dx(d) = dwdx.
Trace();
3351 for (
int d = 0; d <
dim; d++)
3356 Mult(Amat, work2, work1);
3358 d_Winv_dx(d) = work2.
Trace();
3360 d_Winv_dx *= -weight_m;
3362 d_detW_dx += d_Winv_dx;
3424 d_vals.
SetSize(nqp); d_vals = 1.0;
3441 for (
int q = 0; q < nqp; q++)
3467 for (
int i = 0; i < dof; i++)
3469 const double w_shape_i = weight_m * shape(i);
3470 for (
int j = 0; j < dof; j++)
3472 const double w = w_shape_i * shape(j);
3473 for (
int d1 = 0; d1 <
dim; d1++)
3475 for (
int d2 = 0; d2 <
dim; d2++)
3477 elmat(d1*dof + i, d2*dof + j) += w * hess(d1, d2);
3498 Vector shape(dof), adapt_lim_gf_e, adapt_lim_gf_q, adapt_lim_gf0_q(nqp);
3512 grad_phys.
Mult(adapt_lim_gf_e, grad_ptr);
3516 for (
int q = 0; q < nqp; q++)
3520 adapt_lim_gf_grad_e.
MultTranspose(shape, adapt_lim_gf_grad_q);
3521 adapt_lim_gf_grad_q *= 2.0 * (adapt_lim_gf_q(q) - adapt_lim_gf0_q(q));
3534 Vector shape(dof), adapt_lim_gf_e, adapt_lim_gf_q, adapt_lim_gf0_q(nqp);
3548 grad_phys.
Mult(adapt_lim_gf_e, grad_ptr);
3553 Mult(grad_phys, adapt_lim_gf_grad_e, adapt_lim_gf_hess_e);
3560 for (
int q = 0; q < nqp; q++)
3565 adapt_lim_gf_grad_e.
MultTranspose(shape, adapt_lim_gf_grad_q);
3570 for (
int i = 0; i < dof *
dim; i++)
3572 const int idof = i % dof, idim = i / dof;
3573 for (
int j = 0; j <= i; j++)
3575 const int jdof = j % dof, jdim = j / dof;
3576 const double entry =
3577 w * ( 2.0 * adapt_lim_gf_grad_q(idim) * shape(idof) *
3578 adapt_lim_gf_grad_q(jdim) * shape(jdof) +
3579 2.0 * (adapt_lim_gf_q(q) - adapt_lim_gf0_q(q)) *
3580 adapt_lim_gf_hess_q(idim, jdim) * shape(idof) * shape(jdof));
3582 if (i != j) { mat(j, i) += entry; }
3597 for (
int s = 0;
s < dofs.
Size();
s++)
3599 count += ((*surf_fit_marker)[dofs[
s]]) ? 1 : 0;
3601 if (count == 0) {
return; }
3617 grad_phys.Mult(sigma_e, grad_ptr);
3619 Vector shape_x(dof_x), shape_s(dof_s);
3622 surf_fit_grad_s = 0.0;
3624 for (
int s = 0;
s < dof_s;
s++)
3651 for (
int s = 0;
s < dofs.
Size();
s++)
3653 count += ((*surf_fit_marker)[dofs[
s]]) ? 1 : 0;
3655 if (count == 0) {
return; }
3669 grad_phys.Mult(sigma_e, grad_ptr);
3674 Mult(grad_phys, surf_fit_grad_e, surf_fit_hess_e);
3678 Vector shape_x(dof_x), shape_s(dof_s);
3684 for (
int s = 0;
s < dof_s;
s++)
3700 for (
int i = 0; i < dof_x *
dim; i++)
3702 const int idof = i % dof_x, idim = i / dof_x;
3703 for (
int j = 0; j <= i; j++)
3705 const int jdof = j % dof_x, jdim = j / dof_x;
3706 const double entry =
3707 w * ( 2.0 * surf_fit_grad_s(idim) * shape_x(idof) *
3708 surf_fit_grad_s(jdim) * shape_x(jdof) +
3709 2.0 * sigma_e(
s) * surf_fit_hess_s(idim, jdim) *
3710 shape_x(idof) * shape_x(jdof));
3712 if (i != j) { mat(j, i) += entry; }
3720 Vector &elfun,
const int dofidx,
3721 const int dir,
const double e_fx,
3725 int idx = dir*dof+dofidx;
3729 double dfdx = (e_fxph-e_fx)/
dx;
3734 (*(
ElemDer[T.ElementNo]))(idx) = dfdx;
3763 for (
int j = 0; j <
dim; j++)
3765 for (
int i = 0; i < dof; i++)
3795 for (
int q = 0; q < nqp; q++)
3821 for (
int i = 0; i < dof; i++)
3823 for (
int j = 0; j < i+1; j++)
3825 for (
int k1 = 0; k1 <
dim; k1++)
3827 for (
int k2 = 0; k2 <
dim; k2++)
3829 elfunmod(k2*dof+j) +=
dx;
3856 double e_fx = ElemPertLoc(k2*dof+j);
3859 elfunmod(k2*dof+j) -=
dx;
3860 double e_fpx = ElemDerLoc(k1*dof+i);
3862 elmat(k1*dof+i, k2*dof+j) = (e_fpxph - e_fpx) /
dx;
3863 elmat(k2*dof+j, k1*dof+i) = (e_fpxph - e_fpx) /
dx;
3893 for (
int q = 0; q < nqp; q++)
3910 MFEM_VERIFY(cf,
"Dynamic weight works only with a ConstantCoefficient.");
3911 cf->constant *= factor;
3920 MFEM_VERIFY(cf,
"Dynamic weight works only with a ConstantCoefficient.");
3921 return cf->constant;
3950 double &metric_energy,
3952 double &surf_fit_gf_energy)
3963 metric_energy = 0.0;
3965 surf_fit_gf_energy = 0.0;
3966 for (
int i = 0; i <
fes->
GetNE(); i++)
3972 const int dof = fe->
GetDof();
3981 for (
int q = 0; q < nqp; q++)
3986 const double weight = ip.
weight *
Jtr(q).Det();
3993 lim_energy += weight;
4003 for (
int s = 0;
s < dofs.
Size();
s++)
4007 surf_fit_gf_energy += sigma_e(
s) * sigma_e(
s);
4033 dx = std::numeric_limits<float>::max();
4036 double detv_avg_min = std::numeric_limits<float>::max();
4037 for (
int i = 0; i < NE; i++)
4042 for (
int j = 0; j < nsp; j++)
4046 detv_sum += std::fabs(
Jpr.
Det());
4048 double detv_avg = pow(detv_sum/nsp, 1./
dim);
4049 detv_avg_min = std::min(detv_avg, detv_avg_min);
4059 PA.Jtr_needs_update =
true;
4084 MPI_Allreduce(&
dx, &min_jac_all, 1, MPI_DOUBLE, MPI_MIN, pfes->
GetComm());
4097 #ifdef MFEM_USE_GSLIB 4099 if (dynamic_cast<const InterpolatorFP *>(ae))
4101 MFEM_ABORT(
"Using GSLIB-based interpolation with finite differences" 4102 "requires careful consideration. Contact TMOP team.");
4119 #ifdef MFEM_USE_GSLIB 4121 if (dynamic_cast<const InterpolatorFP *>(ae))
4123 MFEM_ABORT(
"Using GSLIB-based interpolation with finite differences" 4124 "requires careful consideration. Contact TMOP team.");
4145 for (
int i = 0; i < NE; i++)
4161 for (
int q = 0; q < nsp; q++)
4170 min_detT = std::min(min_detT, detT);
4196 for (
int i = 0; i < NE; i++)
4215 for (
int q = 0; q < nsp; q++)
4225 double metric_val = 0.0;
4232 max_muT = std::max(max_muT, metric_val);
4244 if (!wcuo) {
return; }
4255 double min_detT_all = min_detT;
4259 MPI_Allreduce(&min_detT, &min_detT_all, 1, MPI_DOUBLE, MPI_MIN,
4263 if (wcuo) { wcuo->
SetMinDetT(min_detT_all); }
4267 double max_muT_all = max_muT;
4271 MPI_Allreduce(&max_muT, &max_muT_all, 1, MPI_DOUBLE, MPI_MAX,
4283 MFEM_VERIFY(
tmopi.Size() > 0,
"No TMOP_Integrators were added.");
4285 tmopi[0]->EnableLimiting(n0, dist, w0, lfunc);
4286 for (
int i = 1; i <
tmopi.Size(); i++) {
tmopi[i]->DisableLimiting(); }
4293 MFEM_VERIFY(
tmopi.Size() > 0,
"No TMOP_Integrators were added.");
4295 tmopi[0]->EnableLimiting(n0, w0, lfunc);
4296 for (
int i = 1; i <
tmopi.Size(); i++) {
tmopi[i]->DisableLimiting(); }
4301 MFEM_VERIFY(
tmopi.Size() > 0,
"No TMOP_Integrators were added.");
4303 tmopi[0]->SetLimitingNodes(n0);
4304 for (
int i = 1; i <
tmopi.Size(); i++) {
tmopi[i]->DisableLimiting(); }
4312 for (
int i = 0; i <
tmopi.Size(); i++)
4314 energy +=
tmopi[i]->GetElementEnergy(el, T, elfun);
4324 MFEM_VERIFY(
tmopi.Size() > 0,
"No TMOP_Integrators were added.");
4326 tmopi[0]->AssembleElementVector(el, T, elfun, elvect);
4327 for (
int i = 1; i <
tmopi.Size(); i++)
4330 tmopi[i]->AssembleElementVector(el, T, elfun, elvect_i);
4340 MFEM_VERIFY(
tmopi.Size() > 0,
"No TMOP_Integrators were added.");
4342 tmopi[0]->AssembleElementGrad(el, T, elfun, elmat);
4343 for (
int i = 1; i <
tmopi.Size(); i++)
4346 tmopi[i]->AssembleElementGrad(el, T, elfun, elmat_i);
4357 for (
int i = 0; i <
tmopi.Size(); i++)
4359 energy +=
tmopi[i]->GetRefinementElementEnergy(el, T, elfun, irule);
4370 for (
int i = 0; i <
tmopi.Size(); i++)
4372 energy +=
tmopi[i]->GetDerefinementElementEnergy(el, T, elfun);
4379 const int cnt =
tmopi.Size();
4380 double total_integral = 0.0;
4381 for (
int i = 0; i < cnt; i++)
4383 tmopi[i]->EnableNormalization(x);
4384 total_integral += 1.0 /
tmopi[i]->metric_normal;
4386 for (
int i = 0; i < cnt; i++)
4388 tmopi[i]->metric_normal = 1.0 / total_integral;
4395 const int cnt =
tmopi.Size();
4396 double total_integral = 0.0;
4397 for (
int i = 0; i < cnt; i++)
4399 tmopi[i]->ParEnableNormalization(x);
4400 total_integral += 1.0 /
tmopi[i]->metric_normal;
4402 for (
int i = 0; i < cnt; i++)
4404 tmopi[i]->metric_normal = 1.0 / total_integral;
4411 for (
int i = 0; i <
tmopi.Size(); i++)
4413 tmopi[i]->AssemblePA(fes);
4420 for (
int i = 0; i <
tmopi.Size(); i++)
4422 tmopi[i]->AssembleGradPA(xe,fes);
4428 for (
int i = 0; i <
tmopi.Size(); i++)
4430 tmopi[i]->AssembleGradDiagonalPA(de);
4436 for (
int i = 0; i <
tmopi.Size(); i++)
4438 tmopi[i]->AddMultPA(xe, ye);
4444 for (
int i = 0; i <
tmopi.Size(); i++)
4446 tmopi[i]->AddMultGradPA(re, ce);
4452 double energy = 0.0;
4453 for (
int i = 0; i <
tmopi.Size(); i++)
4455 energy +=
tmopi[i]->GetLocalStateEnergyPA(xe);
4464 const int NE = mesh.
GetNE();
4472 for (
int i = 0; i < NE; i++)
4489 for (
int j = 0; j < nsp; j++)
4500 metric_gf(gf_dofs[j]) = metric.
EvalW(T);
Abstract class for all finite elements.
void ComputeUntangleMetricQuantiles(const Vector &x, const FiniteElementSpace &fes)
InvariantsEvaluator2D< double > ie
virtual void ComputeElementTargets(int e_id, const FiniteElement &fe, const IntegrationRule &ir, const Vector &elfun, DenseTensor &Jtr) const
Given an element and quadrature rule, computes ref->target transformation Jacobians for each quadratu...
virtual double GetDerefinementElementEnergy(const FiniteElement &el, ElementTransformation &T, const Vector &elfun)
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
virtual double EvalWMatrixForm(const DenseMatrix &Jpt) const
Evaluates the metric in matrix form (opposed to invariant form). Used for validating the invariant ev...
virtual void SetParDiscreteTargetSize(const ParGridFunction &tspec_)
virtual void ComputeElementTargets(int e_id, const FiniteElement &fe, const IntegrationRule &ir, const Vector &elfun, DenseTensor &Jtr) const
Given an element and quadrature rule, computes ref->target transformation Jacobians for each quadratu...
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
InvariantsEvaluator3D< double > ie
virtual void SetInitialField(const Vector &init_nodes, const Vector &init_field)=0
int GetNPoints() const
Returns the number of the points in the integration rule.
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const =0
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
InvariantsEvaluator3D< double > ie
void AssembleElemVecSurfFit(const FiniteElement &el_x, IsoparametricTransformation &Tpr, DenseMatrix &mat)
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
Class for an integration rule - an Array of IntegrationPoint.
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
VectorCoefficient * vector_tspec
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
void EnableFiniteDifferences(const GridFunction &x)
Enables FD-based approximation and computes dx.
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
const TargetType target_type
Class for grid function - Vector with associated FE space.
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const double weight, DenseMatrix &A) const
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
void ComputeAvgMetrics(const GridFunction &nodes, const TargetConstructor &tc, Vector &averages) const
void UseExternalData(double *ext_data, int i, int j, int k)
virtual void Update(bool want_transform=true)
Reflect changes in the mesh: update number of DOFs, etc. Also, calculate GridFunction transformation ...
InvariantsEvaluator3D< double > ie
InvariantsEvaluator3D< double > ie
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
void Assemble_ddI2b(scalar_t w, scalar_t *A)
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
double GetSurfaceFittingWeight()
Get the surface fitting weight.
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
Base class for vector Coefficients that optionally depend on time and space.
A coefficient that is constant across space and time.
virtual double EvalW(const DenseMatrix &Jpt) const
Evaluate the strain energy density function, W = W(Jpt), by using the 2D or 3D matrix invariants...