Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "UELThread.h"
11 : #include "NonlinearSystemBase.h"
12 :
13 908 : UELThread::UELThread(FEProblemBase & fe_problem, AbaqusUserElement & uel_uo)
14 : : ThreadedElementLoop<ConstElemRange>(fe_problem),
15 908 : _sys(fe_problem.currentNonlinearSystem()),
16 908 : _aux_sys(&_fe_problem.getAuxiliarySystem()),
17 908 : _variables(uel_uo.getVariables()),
18 908 : _aux_variables(uel_uo.getAuxVariables()),
19 908 : _uel_uo(uel_uo),
20 908 : _uel(uel_uo.getPlugin()),
21 908 : _statev_copy(uel_uo._nstatev)
22 : {
23 : // general solve step
24 908 : _lflags[3] = 0;
25 :
26 : // newton based solution (should be 1 when using a predictor or maybe at the start of the
27 : // iteration)
28 908 : _lflags[4] = 0;
29 908 : }
30 :
31 : // Splitting Constructor
32 0 : UELThread::UELThread(UELThread & x, Threads::split split)
33 : : ThreadedElementLoop<ConstElemRange>(x, split),
34 0 : _sys(x._sys),
35 0 : _aux_sys(x._aux_sys),
36 0 : _variables(x._variables),
37 0 : _aux_variables(x._aux_variables),
38 0 : _uel_uo(x._uel_uo),
39 0 : _uel(x._uel),
40 0 : _lflags(x._lflags),
41 0 : _statev_copy(_uel_uo._nstatev)
42 : {
43 0 : }
44 :
45 : void
46 908 : UELThread::subdomainChanged()
47 : {
48 908 : }
49 :
50 : /**
51 : * Fortran array memory layout: (1,1), (2,1) (3,1) (1,2) (2,2) (3,2) (1,3) (2,3) (3,3)
52 : * C++ array memory layout: [0][0], [0][1], [0][2], [1][0], [1][1], [1][2], [2][0], [2][1],
53 : * [2][2]
54 : */
55 : void
56 240216 : UELThread::onElement(const Elem * elem)
57 : {
58 : // get all dof indices for the coupled variables on this element
59 240216 : int nnode = elem->n_nodes();
60 240216 : const auto nvar = _variables.size();
61 240216 : _all_dof_indices.resize(nnode * nvar);
62 :
63 : // ** System variables **
64 : // prepare DOF indices in the correct order
65 720864 : for (const auto i : index_range(_variables))
66 : {
67 480648 : _variables[i]->getDofIndices(elem, _var_dof_indices);
68 :
69 : // sanity checks
70 480648 : if (_variables[i]->scalingFactor() != 1)
71 0 : mooseError("Scaling factors other than unity are not yet supported");
72 480648 : if (static_cast<int>(_var_dof_indices.size()) != nnode)
73 0 : mooseError("All coupled variables must be full order lagrangian");
74 :
75 1925832 : for (const auto j : make_range(nnode))
76 1445184 : _all_dof_indices[j * nvar + i] = _var_dof_indices[j];
77 : }
78 :
79 240216 : int ndofel = _all_dof_indices.size();
80 :
81 : // Get solution values
82 240216 : _all_dof_values.resize(ndofel);
83 :
84 240216 : _sys.currentSolution()->get(_all_dof_indices, _all_dof_increments);
85 240216 : _all_dof_increments.resize(ndofel);
86 240216 : _sys.solutionOld().get(_all_dof_indices, _all_dof_values);
87 :
88 : mooseAssert(_all_dof_values.size() == _all_dof_increments.size(), "Inconsistent solution size.");
89 1685400 : for (const auto i : index_range(_all_dof_values))
90 1445184 : _all_dof_increments[i] -= _all_dof_values[i];
91 :
92 240216 : _all_udot_dof_values.resize(ndofel);
93 240216 : _all_udotdot_dof_values.resize(ndofel);
94 :
95 : // Get u_dot and u_dotdot solution values (required for future expansion of the interface)
96 : if (false)
97 : {
98 : _all_udot_dof_values.resize(ndofel);
99 : _sys.solutionUDot()->get(_all_dof_indices, _all_udot_dof_values);
100 : _all_udotdot_dof_values.resize(ndofel);
101 : _sys.solutionUDot()->get(_all_dof_indices, _all_udotdot_dof_values);
102 : }
103 :
104 : // Prepare external fields
105 240216 : const auto nvar_aux = _aux_variables.size();
106 240216 : _all_aux_var_dof_indices.resize(nnode * nvar_aux);
107 240216 : _all_aux_var_dof_increments.resize(nnode * nvar_aux);
108 :
109 432648 : for (const auto i : index_range(_aux_variables))
110 : {
111 192432 : _aux_variables[i]->getDofIndices(elem, _aux_var_dof_indices);
112 :
113 192432 : if (static_cast<int>(_aux_var_dof_indices.size()) != nnode)
114 0 : mooseError("All auxiliary variables must be full order Lagrangian");
115 :
116 771888 : for (const auto j : make_range(nnode))
117 579456 : _all_aux_var_dof_indices[j * nvar_aux + i] = _aux_var_dof_indices[j];
118 : }
119 :
120 240216 : _all_aux_var_dof_values.resize(nnode * nvar_aux);
121 240216 : _aux_var_values_to_uel.resize(nnode * nvar_aux * 2); // Value _and_ increment
122 :
123 240216 : _aux_sys->currentSolution()->get(_all_aux_var_dof_indices, _all_aux_var_dof_increments);
124 240216 : _all_aux_var_dof_increments.resize(nnode * nvar_aux);
125 240216 : _aux_sys->solutionOld().get(_all_aux_var_dof_indices, _all_aux_var_dof_values);
126 :
127 : // First, one external field and increment; then, second external field and increment, etc.
128 819672 : for (const auto i : index_range(_all_aux_var_dof_values))
129 579456 : _all_aux_var_dof_increments[i] -= _all_aux_var_dof_values[i];
130 :
131 819672 : for (const auto i : index_range(_all_aux_var_dof_values))
132 : {
133 579456 : _aux_var_values_to_uel[i * 2] = _all_aux_var_dof_values[i];
134 579456 : _aux_var_values_to_uel[i * 2 + 1] = _all_aux_var_dof_increments[i];
135 : }
136 : // End of prepare external fields.
137 :
138 240216 : const bool do_residual = _sys.hasVector(_sys.residualVectorTag());
139 240216 : const bool do_jacobian = _sys.hasMatrix(_sys.systemMatrixTag());
140 :
141 : // set flags
142 240216 : if (do_residual && do_jacobian)
143 0 : _lflags[2] = 0;
144 240216 : else if (!do_residual && do_jacobian)
145 60054 : _lflags[2] = 4;
146 180162 : else if (do_residual && !do_jacobian)
147 180162 : _lflags[2] = 5;
148 : else
149 0 : _lflags[2] = 0;
150 :
151 : // copy coords
152 240216 : int dim = _uel_uo._dim;
153 240216 : _coords.resize(nnode * dim);
154 961944 : for (const auto i : make_range(nnode))
155 2166912 : for (const auto j : make_range(dim))
156 1445184 : _coords[j + dim * i] = elem->node_ref(i)(j);
157 :
158 : // prepare residual and jacobian storage
159 240216 : _local_re.resize(ndofel);
160 240216 : _local_ke.resize(ndofel, ndofel);
161 :
162 240216 : int nrhs = 1; // : RHS should contain the residual vector
163 240216 : int jtype = _uel_uo._jtype; // type of user element
164 :
165 240216 : Real dt = _fe_problem.dt();
166 240216 : Real time = _fe_problem.time();
167 240216 : std::vector<Real> times{time - dt, time - dt}; // first entry should be the step time (TODO)
168 :
169 : std::array<Real, 8> energy;
170 240216 : int jelem = elem->id() + 1; // User-assigned element number
171 : Real pnewdt;
172 :
173 240216 : int npredf = nvar_aux; // Number of external fields.
174 :
175 : // dummy vars
176 240216 : Real rdummy = 0;
177 240216 : int idummy = 0;
178 :
179 : // stateful data
180 240216 : if (_uel_uo._nstatev)
181 : {
182 240216 : const auto & statev_old = _uel_uo._statev[_uel_uo._statev_index_old][elem->id()];
183 240216 : std::copy(statev_old.begin(), statev_old.end(), _statev_copy.begin());
184 : }
185 :
186 : // call the plugin
187 240216 : _uel(_local_re.get_values().data(),
188 : _local_ke.get_values().data(),
189 : _statev_copy.data(),
190 : energy.data(),
191 : &ndofel,
192 : &nrhs,
193 : &_uel_uo._nstatev,
194 : _uel_uo._props.data(),
195 240216 : &_uel_uo._nprops,
196 : _coords.data(),
197 : &dim,
198 : &nnode,
199 : _all_dof_values.data() /* U[] */,
200 : _all_dof_increments.data() /* DU[] */,
201 : _all_udot_dof_values.data() /* V[] */,
202 : _all_udotdot_dof_values.data() /* A[] */,
203 : &jtype,
204 : times.data(),
205 : &dt,
206 : &idummy /* KSTEP */,
207 : &idummy /* KINC */,
208 : &jelem,
209 : nullptr /* PARAMS[] */,
210 : &idummy /* NDLOAD */,
211 : nullptr /* JDLTYP[] */,
212 : nullptr /* ADLMAG[] */,
213 : _aux_var_values_to_uel.data() /* PREDEF[] */,
214 : &npredf /* NPREDF */,
215 : nullptr /* LFLAGS[] */,
216 : &ndofel /* MLVARX */,
217 : nullptr /* DDLMAG[] */,
218 : &idummy /* MDLOAD */,
219 : &pnewdt,
220 : nullptr /* JPROPS[] */,
221 : &idummy /* NJPROP */,
222 : &rdummy /* PERIOD */
223 : );
224 :
225 240216 : if (_uel_uo._nstatev)
226 : {
227 240216 : auto & statev_current = _uel_uo._statev[_uel_uo._statev_index_current][elem->id()];
228 240216 : std::copy(_statev_copy.begin(), _statev_copy.end(), statev_current.begin());
229 : }
230 :
231 : // write to the residual vector
232 : // sign of 'residuals' has been tested with external loading and matches that of moose-umat
233 : // setups.
234 240216 : if (do_residual)
235 180162 : _uel_uo.addResiduals(
236 180162 : _fe_problem.assembly(_tid, _sys.number()), _local_re, _all_dof_indices, -1.0);
237 :
238 : // write to the Jacobian (hope we don't have to transpose...)
239 240216 : if (do_jacobian)
240 60054 : _uel_uo.addJacobian(_fe_problem.assembly(_tid, _sys.number()),
241 : _local_ke,
242 : _all_dof_indices,
243 : _all_dof_indices,
244 : -1.0);
245 240216 : }
|
