Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 1424 : UELThread::UELThread(FEProblemBase & fe_problem, AbaqusUserElement & uel_uo)
14 : : ThreadedElementLoop<ConstElemRange>(fe_problem),
15 1424 : _sys(fe_problem.currentNonlinearSystem()),
16 1424 : _aux_sys(&_fe_problem.getAuxiliarySystem()),
17 1424 : _variables(uel_uo.getVariables()),
18 1424 : _aux_variables(uel_uo.getAuxVariables()),
19 1424 : _uel_uo(uel_uo),
20 1424 : _uel(uel_uo.getPlugin()),
21 1424 : _statev_copy(uel_uo._nstatev)
22 : {
23 : // general solve step
24 1424 : _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 1424 : _lflags[4] = 0;
29 1424 : }
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 1424 : UELThread::subdomainChanged()
47 : {
48 1424 : }
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 361998 : UELThread::onElement(const Elem * elem)
57 : {
58 : // get all dof indices for the coupled variables on this element
59 361998 : int nnode = elem->n_nodes();
60 361998 : const auto nvar = _variables.size();
61 361998 : _all_dof_indices.resize(nnode * nvar);
62 :
63 : // ** System variables **
64 : // prepare DOF indices in the correct order
65 1087992 : for (const auto i : index_range(_variables))
66 : {
67 725994 : _variables[i]->getDofIndices(elem, _var_dof_indices);
68 :
69 : // sanity checks
70 725994 : if (_variables[i]->scalingFactor() != 1)
71 0 : mooseError("Scaling factors other than unity are not yet supported");
72 725994 : if (static_cast<int>(_var_dof_indices.size()) != nnode)
73 0 : mooseError("All coupled variables must be full order lagrangian");
74 :
75 2933946 : for (const auto j : make_range(nnode))
76 2207952 : _all_dof_indices[j * nvar + i] = _var_dof_indices[j];
77 : }
78 :
79 361998 : int ndofel = _all_dof_indices.size();
80 :
81 : // Get solution values
82 361998 : _all_dof_values.resize(ndofel);
83 :
84 361998 : _sys.currentSolution()->get(_all_dof_indices, _all_dof_increments);
85 361998 : _all_dof_increments.resize(ndofel);
86 361998 : _sys.solutionOld().get(_all_dof_indices, _all_dof_values);
87 :
88 : mooseAssert(_all_dof_values.size() == _all_dof_increments.size(), "Inconsistent solution size.");
89 2569950 : for (const auto i : index_range(_all_dof_values))
90 2207952 : _all_dof_increments[i] -= _all_dof_values[i];
91 :
92 361998 : _all_udot_dof_values.resize(ndofel);
93 361998 : _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 361998 : const auto nvar_aux = _aux_variables.size();
106 361998 : _all_aux_var_dof_indices.resize(nnode * nvar_aux);
107 361998 : _all_aux_var_dof_increments.resize(nnode * nvar_aux);
108 :
109 653994 : for (const auto i : index_range(_aux_variables))
110 : {
111 291996 : _aux_variables[i]->getDofIndices(elem, _aux_var_dof_indices);
112 :
113 291996 : if (static_cast<int>(_aux_var_dof_indices.size()) != nnode)
114 0 : mooseError("All auxiliary variables must be full order Lagrangian");
115 :
116 1187964 : for (const auto j : make_range(nnode))
117 895968 : _all_aux_var_dof_indices[j * nvar_aux + i] = _aux_var_dof_indices[j];
118 : }
119 :
120 361998 : _all_aux_var_dof_values.resize(nnode * nvar_aux);
121 361998 : _aux_var_values_to_uel.resize(nnode * nvar_aux * 2); // Value _and_ increment
122 :
123 361998 : _aux_sys->currentSolution()->get(_all_aux_var_dof_indices, _all_aux_var_dof_increments);
124 361998 : _all_aux_var_dof_increments.resize(nnode * nvar_aux);
125 361998 : _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 1257966 : for (const auto i : index_range(_all_aux_var_dof_values))
129 895968 : _all_aux_var_dof_increments[i] -= _all_aux_var_dof_values[i];
130 :
131 1257966 : for (const auto i : index_range(_all_aux_var_dof_values))
132 : {
133 895968 : _aux_var_values_to_uel[i * 2] = _all_aux_var_dof_values[i];
134 895968 : _aux_var_values_to_uel[i * 2 + 1] = _all_aux_var_dof_increments[i];
135 : }
136 : // End of prepare external fields.
137 :
138 361998 : const bool do_residual = _sys.hasVector(_sys.residualVectorTag());
139 361998 : const bool do_jacobian = _sys.hasMatrix(_sys.systemMatrixTag());
140 :
141 : // set flags
142 361998 : if (do_residual && do_jacobian)
143 0 : _lflags[2] = 0;
144 361998 : else if (!do_residual && do_jacobian)
145 120108 : _lflags[2] = 4;
146 241890 : else if (do_residual && !do_jacobian)
147 241890 : _lflags[2] = 5;
148 : else
149 0 : _lflags[2] = 0;
150 :
151 : // copy coords
152 361998 : int dim = _uel_uo._dim;
153 361998 : _coords.resize(nnode * dim);
154 1457982 : for (const auto i : make_range(nnode))
155 3303936 : for (const auto j : make_range(dim))
156 2207952 : _coords[j + dim * i] = elem->node_ref(i)(j);
157 :
158 : // prepare residual and jacobian storage
159 361998 : _local_re.resize(ndofel);
160 361998 : _local_ke.resize(ndofel, ndofel);
161 :
162 361998 : int nrhs = 1; // : RHS should contain the residual vector
163 361998 : int jtype = _uel_uo._jtype; // type of user element
164 :
165 361998 : Real dt = _fe_problem.dt();
166 361998 : Real time = _fe_problem.time();
167 361998 : std::vector<Real> times{time - dt, time - dt}; // first entry should be the step time (TODO)
168 :
169 : std::array<Real, 8> energy;
170 361998 : int jelem = elem->id() + 1; // User-assigned element number
171 : Real pnewdt;
172 :
173 361998 : int npredf = nvar_aux; // Number of external fields.
174 :
175 : // dummy vars
176 361998 : Real rdummy = 0;
177 361998 : int idummy = 0;
178 :
179 : // stateful data
180 361998 : if (_uel_uo._nstatev)
181 : {
182 361998 : const auto & statev_old = _uel_uo._statev[_uel_uo._statev_index_old][elem->id()];
183 361998 : std::copy(statev_old.begin(), statev_old.end(), _statev_copy.begin());
184 : }
185 :
186 : // call the plugin
187 361998 : _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 361998 : &_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 361998 : if (_uel_uo._nstatev)
226 : {
227 361998 : auto & statev_current = _uel_uo._statev[_uel_uo._statev_index_current][elem->id()];
228 361998 : 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 361998 : if (do_residual)
235 241890 : _uel_uo.addResiduals(
236 241890 : _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 361998 : if (do_jacobian)
240 120108 : _uel_uo.addJacobian(_fe_problem.assembly(_tid, _sys.number()),
241 : _local_ke,
242 : _all_dof_indices,
243 : _all_dof_indices,
244 : -1.0);
245 361998 : }
|
