1!
2!     CalculiX - A 3-dimensional finite element program
3!     Copyright (C) 1998-2021 Guido Dhondt
4!
5!     This program is free software; you can redistribute it and/or
6!     modify it under the terms of the GNU General Public License as
7!     published by the Free Software Foundation(version 2);
8!
9!
10!     This program is distributed in the hope that it will be useful,
11!     but WITHOUT ANY WARRANTY; without even the implied warranty of
12!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13!     GNU General Public License for more details.
14!
15!     You should have received a copy of the GNU General Public License
16!     along with this program; if not, write to the Free Software
17!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18!
19!     construction of the B matrix
20!
21      subroutine networkmpc_lhs(i,ipompc,nodempc,coefmpc,
22     &  labmpc,v,nactdog,ac,j,mi,nteq,ipkon,kon,lakon,iponoel,
23     &  inoel,ielprop,prop,ielmat,
24     &  shcon,nshcon,rhcon,nrhcon,ntmat_,cocon,ncocon)
25!
26!     user defined network mpc: calculation of the left hand
27!     side
28!
29!     INPUT:
30!
31!     i                  MPC number
32!     ipompc(1..nmpc))   ipompc(i) points to the first term of
33!                        MPC i in field nodempc
34!     nodempc(1,*)       node number of a MPC term
35!     nodempc(2,*)       coordinate direction of a MPC term
36!     nodempc(3,*)       if not 0: points towards the next term
37!                                  of the MPC in field nodempc
38!                        if 0: MPC definition is finished
39!     coefmpc(*)         coefficient of a MPC term
40!     labmpc(*)          label of the MPC. For user-defined
41!                        network MPC's it starts with NETWORK;
42!                        the remaining 13 characters can be used
43!                        to distinguish between different kinds of
44!                        network user MPC's
45!     v(0..mi(2),1..nk)  actual solution field in all nodes
46!                        0: total temperature
47!                        1: mass flow
48!                        2: total pressure
49!     nactdog(j,i)       determines the network equation corresponding
50!                        to degree of freedom j in node i;
51!                        if zero the degree of freedom is not active
52!     j                  network equation corresponding to the
53!                        present MPC (i.e. MPC i)
54!     mi(*)              field with global information; mi(2) is the
55!                        highest variable number
56!     nteq               number of network equations
57!     ipkon(i)           points to the location in field kon preceding
58!                        the topology of element i
59!     kon(*)             contains the topology of all elements. The
60!                        topology of element i starts at kon(ipkon(i)+1)
61!                        and continues until all nodes are covered. The
62!                        number of nodes depends on the element label
63!     lakon(i)           contains the label of element i
64!     iponoel(i)         the network elements to which node i belongs
65!                        are stored in inoel(1,iponoel(i)),
66!                        inoel(1,inoel(2,iponoel(i)))...... until
67!                        inoel(2,inoel(2,inoel(2......)=0
68!     inoel(1..2,*)      field containing the network elements
69!     ielprop(i)         points to the location in field prop preceding
70!                        the properties of element i
71!     prop(*)            contains the properties of all network elements. The
72!                        properties of element i start at prop(ielprop(i)+1)
73!                        and continues until all properties are covered. The
74!                        appropriate amount of properties depends on the
75!                        element label. The kind of properties, their
76!                        number and their order corresponds
77!                        to the description in the user's manual,
78!                        cf. the sections "Fluid Section Types"
79!     ielmat(j,i)        contains the material number for element i
80!                        and layer j
81!     shcon(0,j,i)       temperature at temperature point j of material i
82!     shcon(1,j,i)       specific heat at constant pressure at the
83!                        temperature point j of material i
84!     shcon(2,j,i)       dynamic viscosity at the temperature point j of
85!                        material i
86!     shcon(3,1,i)       specific gas constant of material i
87!     nshcon(i)          number of temperature data points for the specific
88!                        heat of material i
89!     rhcon(0,j,i)       temperature at density temperature point j of
90!                        material i
91!     rhcon(1,j,i)       density at the density temperature point j of
92!                        material i
93!     nrhcon(i)          number of temperature data points for the density
94!                        of material i
95!     ntmat_             maximum number of temperature data points for
96!                        any material property for any material
97!     ncocon(1,i)        number of conductivity constants for material i
98!     ncocon(2,i)        number of temperature data points for the
99!                        conductivity coefficients of material i
100!     cocon(0,j,i)       temperature at conductivity temperature point
101!                        j of material i
102!     cocon(k,j,i)       conductivity coefficient k at conductivity
103!                        temperature point j of material i
104!
105!     OUTPUT:
106!
107!     ac(*)              left hand side of the system of network
108!                        equations; this routines should return the
109!                        derivative of the network MPC at stake w.r.t.
110!                        all active degrees of freedom occurring in the
111!                        MPC and store them in row j of matrix ac.
112!
113      implicit none
114!
115      character*8 lakon(*)
116      character*20 labmpc(*)
117!
118      integer mi(*),i,ipompc(*),nodempc(3,*),j,index,node,idir,
119     &  nactdog(0:3,*),nteq,ipkon(*),kon(*),
120     &  iponoel(*),inoel(2,*),ielprop(*),ielmat(mi(3),*),nshcon(*),
121     &  nrhcon(*),ncocon(2,*),ntmat_
122!
123      real*8 coefmpc(*),v(0:mi(2),*),ac(nteq,*),
124     &  prop(*),shcon(0:3,ntmat_,*),rhcon(0:1,ntmat_,*),
125     &  cocon(0:6,ntmat_,*)
126!
127!
128!
129      if(labmpc(i)(8:16).eq.'QUADRATIC') then
130!
131!        example equation of the form
132!        f:=a*v(idir1,node1)+b*v(idir2,node2)**2=0
133!
134!        a,idir1,node1,b,idir2,node2 are given in the input deck
135!        using the *NETWORK MPC keyword
136!        to be calculated: a*df/d(v(idir1,node1))
137!                          b*df/d(v(idir2,node2))
138!
139         index=ipompc(i)
140         node=nodempc(1,index)
141         idir=nodempc(2,index)
142         ac(j,nactdog(idir,node))=coefmpc(index)
143!
144         index=nodempc(3,index)
145         node=nodempc(1,index)
146         idir=nodempc(2,index)
147!
148!        if nactdog(idir,node) is zero the degree of freedom is
149!        not active
150!
151         if(nactdog(idir,node).ne.0) then
152            ac(j,nactdog(idir,node))=2.d0*coefmpc(index)*v(idir,node)
153         endif
154      else
155         write(*,*) '*ERROR in networkmpc_lhs:'
156         write(*,*) '       unknown MPC: ',labmpc(i)
157         call exit(201)
158      endif
159!
160      return
161      end
162