diff --git a/docs/theory/unsrtnat.bst b/docs/theory/apalike.bst
similarity index 60%
copy from docs/theory/unsrtnat.bst
copy to docs/theory/apalike.bst
index 1c7eddf..65d2585 100644
--- a/docs/theory/unsrtnat.bst
+++ b/docs/theory/apalike.bst
@@ -1,1343 +1,1105 @@
-%% File: `unsrtnat.bst'
-%% A modification of `unsrt.bst' for use with natbib package
-%%
-%% Copyright 1993-2007 Patrick W Daly
-%% Max-Planck-Institut f\"ur Sonnensystemforschung
-%% Max-Planck-Str. 2
-%% D-37191 Katlenburg-Lindau
-%% Germany
-%% E-mail: daly@mps.mpg.de
-%%
-%% This program can be redistributed and/or modified under the terms
-%% of the LaTeX Project Public License Distributed from CTAN
-%% archives in directory macros/latex/base/lppl.txt; either
-%% version 1 of the License, or any later version.
-%%
- % Version and source file information:
- % \ProvidesFile{natbst.mbs}[2007/11/26 1.93 (PWD)]
- %
- % BibTeX `plainnat' family
- % version 0.99b for BibTeX versions 0.99a or later,
- % for LaTeX versions 2.09 and 2e.
- %
- % For use with the `natbib.sty' package; emulates the corresponding
- % member of the `plain' family, but with author-year citations.
- %
- % With version 6.0 of `natbib.sty', it may also be used for numerical
- % citations, while retaining the commands \citeauthor, \citefullauthor,
- % and \citeyear to print the corresponding information.
- %
- % For version 7.0 of `natbib.sty', the KEY field replaces missing
- % authors/editors, and the date is left blank in \bibitem.
- %
- % Includes field EID for the sequence/citation number of electronic journals
- % which is used instead of page numbers.
- %
- % Includes fields ISBN and ISSN.
- %
- % Includes field URL for Internet addresses.
- %
- % Includes field DOI for Digital Object Idenfifiers.
- %
- % Works best with the url.sty package of Donald Arseneau.
- %
- % Works with identical authors and year are further sorted by
- % title text, as in the standard plain.bst etc.
- %
+% BibTeX `apalike' bibliography style (version 0.99a, 8-Dec-10), adapted from
+% the `alpha' style, version 0.99a; for BibTeX version 0.99a.
+%
+% Copyright (C) 1988, 2010 Oren Patashnik.
+% Unlimited copying and redistribution of this file are permitted as long as
+% it is unmodified. Modifications (and redistribution of modified versions)
+% are also permitted, but only if the resulting file is renamed.
+%
+% Differences between this style and `alpha' are generally heralded by a `%'.
+% The file btxbst.doc has the documentation for alpha.bst.
+%
+% This style should be used with the `apalike' LaTeX style (apalike.sty).
+% \cite's come out like "(Jones, 1986)" in the text but there are no labels
+% in the bibliography, and something like "(1986)" comes out immediately
+% after the author. Author (and editor) names appear as last name, comma,
+% initials. A `year' field is required for every entry, and so is either
+% an author (or in some cases, an editor) field or a key field.
+%
+% Editorial note:
+% Many journals require a style like `apalike', but I strongly, strongly,
+% strongly recommend that you not use it if you have a choice---use something
+% like `plain' instead. Mary-Claire van Leunen (A Handbook for Scholars,
+% Knopf, 1979) argues convincingly that a style like `plain' encourages better
+% writing than one like `apalike'. Furthermore the strongest arguments for
+% using an author-date style like `apalike'---that it's "the most practical"
+% (The Chicago Manual of Style, University of Chicago Press, thirteenth
+% edition, 1982, pages 400--401)---fall flat on their face with the new
+% computer-typesetting technology. For instance page 401 anachronistically
+% states "The chief disadvantage of [a style like `plain'] is that additions
+% or deletions cannot be made after the manuscript is typed without changing
+% numbers in both text references and list." LaTeX sidesteps the disadvantage.
+%
+% History:
+% 15-sep-86 (OP) Original version by Oren Patashnik, ideas from Susan King.
+% 10-nov-86 (OP) Truncated the sort.key$ string to the correct length
+% in bib.sort.order to eliminate error message.
+% 24-jan-88 (OP) Updated for BibTeX version 0.99a, from alpha.bst 0.99a;
+% apalike now sorts by author, then year, then title;
+% THIS `apalike' VERSION DOES NOT WORK WITH BIBTEX 0.98i.
+% 8-dec-10 (OP) Still version 0.99a, as the code itself was unchanged;
+% this release clarified the license.
+
ENTRY
{ address
author
booktitle
chapter
- doi
- eprint
- eid
edition
editor
howpublished
institution
- isbn
- issn
journal
key
- month
+% month not used in apalike
note
number
organization
pages
publisher
school
series
title
type
- url
volume
year
}
{}
- { label extra.label sort.label short.list }
+ { label extra.label sort.label }
INTEGERS { output.state before.all mid.sentence after.sentence after.block }
FUNCTION {init.state.consts}
{ #0 'before.all :=
#1 'mid.sentence :=
#2 'after.sentence :=
#3 'after.block :=
}
STRINGS { s t }
FUNCTION {output.nonnull}
{ 's :=
output.state mid.sentence =
{ ", " * write$ }
{ output.state after.block =
{ add.period$ write$
newline$
"\newblock " write$
}
{ output.state before.all =
'write$
{ add.period$ " " * write$ }
if$
}
if$
mid.sentence 'output.state :=
}
if$
s
}
FUNCTION {output}
{ duplicate$ empty$
'pop$
'output.nonnull
if$
}
FUNCTION {output.check}
{ 't :=
duplicate$ empty$
{ pop$ "empty " t * " in " * cite$ * warning$ }
'output.nonnull
if$
}
+% apalike needs this function because
+% the year has special punctuation;
+% apalike ignores the month
+FUNCTION {output.year.check}
+{ year empty$
+ { "empty year in " cite$ * warning$ }
+ { write$
+ " (" year * extra.label * ")" *
+ mid.sentence 'output.state :=
+ }
+ if$
+}
+
+FUNCTION {output.bibitem}
+{ newline$
+ "\bibitem[" write$
+ label write$
+ "]{" write$
+ cite$ write$
+ "}" write$
+ newline$
+ ""
+ before.all 'output.state :=
+}
+
FUNCTION {fin.entry}
{ add.period$
write$
newline$
}
FUNCTION {new.block}
{ output.state before.all =
'skip$
{ after.block 'output.state := }
if$
}
FUNCTION {new.sentence}
{ output.state after.block =
'skip$
{ output.state before.all =
'skip$
{ after.sentence 'output.state := }
if$
}
if$
}
FUNCTION {not}
{ { #0 }
{ #1 }
if$
}
FUNCTION {and}
{ 'skip$
{ pop$ #0 }
if$
}
FUNCTION {or}
{ { pop$ #1 }
'skip$
if$
}
-FUNCTION {new.block.checka}
-{ empty$
- 'skip$
- 'new.block
- if$
-}
-
FUNCTION {new.block.checkb}
{ empty$
swap$ empty$
and
'skip$
'new.block
if$
}
-FUNCTION {new.sentence.checka}
-{ empty$
- 'skip$
- 'new.sentence
- if$
-}
-
-FUNCTION {new.sentence.checkb}
-{ empty$
- swap$ empty$
- and
- 'skip$
- 'new.sentence
- if$
-}
-
FUNCTION {field.or.null}
{ duplicate$ empty$
{ pop$ "" }
'skip$
if$
}
FUNCTION {emphasize}
{ duplicate$ empty$
{ pop$ "" }
- { "\emph{" swap$ * "}" * }
+ { "{\em " swap$ * "}" * }
if$
}
INTEGERS { nameptr namesleft numnames }
FUNCTION {format.names}
{ 's :=
#1 'nameptr :=
s num.names$ 'numnames :=
numnames 'namesleft :=
{ namesleft #0 > }
- { s nameptr "{ff~}{vv~}{ll}{, jj}" format.name$ 't :=
+ { s nameptr "{vv~}{ll}{, jj}{, ff~}" format.name$ 't := % last name first
nameptr #1 >
{ namesleft #1 >
{ ", " * t * }
{ numnames #2 >
{ "," * }
'skip$
if$
t "others" =
{ " et~al." * }
{ " and " * t * }
if$
}
if$
}
't
if$
nameptr #1 + 'nameptr :=
namesleft #1 - 'namesleft :=
}
while$
}
-FUNCTION {format.key}
-{ empty$
- { key field.or.null }
+FUNCTION {format.authors}
+{ author empty$
{ "" }
+ { author format.names }
if$
}
-FUNCTION {format.authors}
-{ author empty$
+FUNCTION {format.key} % this function is just for apalike
+{ empty$
+ { key field.or.null }
{ "" }
- { author format.names }
if$
}
FUNCTION {format.editors}
{ editor empty$
{ "" }
{ editor format.names
editor num.names$ #1 >
{ ", editors" * }
{ ", editor" * }
if$
}
if$
}
-FUNCTION {format.isbn}
-{ isbn empty$
- { "" }
- { new.block "ISBN " isbn * }
- if$
-}
-
-FUNCTION {format.issn}
-{ issn empty$
- { "" }
- { new.block "ISSN " issn * }
- if$
-}
-
-FUNCTION {format.url}
-{ url empty$
- { "" }
- { new.block "URL \url{" url * "}" * }
- if$
-}
-
-FUNCTION {format.doi}
-{ doi empty$
- { "" }
- { new.block "\doi{" doi * "}" * }
- if$
-}
-
-FUNCTION {format.eprint}
-{ eprint empty$
- { "" }
- { new.block "\eprint{" eprint * "}" * }
- if$
-}
-
FUNCTION {format.title}
{ title empty$
{ "" }
{ title "t" change.case$ }
if$
}
-FUNCTION {format.full.names}
-{'s :=
- #1 'nameptr :=
- s num.names$ 'numnames :=
- numnames 'namesleft :=
- { namesleft #0 > }
- { s nameptr
- "{vv~}{ll}" format.name$ 't :=
- nameptr #1 >
- {
- namesleft #1 >
- { ", " * t * }
- {
- numnames #2 >
- { "," * }
- 'skip$
- if$
- t "others" =
- { " et~al." * }
- { " and " * t * }
- if$
- }
- if$
- }
- 't
- if$
- nameptr #1 + 'nameptr :=
- namesleft #1 - 'namesleft :=
- }
- while$
-}
-
-FUNCTION {author.editor.full}
-{ author empty$
- { editor empty$
- { "" }
- { editor format.full.names }
- if$
- }
- { author format.full.names }
- if$
-}
-
-FUNCTION {author.full}
-{ author empty$
- { "" }
- { author format.full.names }
- if$
-}
-
-FUNCTION {editor.full}
-{ editor empty$
- { "" }
- { editor format.full.names }
- if$
-}
-
-FUNCTION {make.full.names}
-{ type$ "book" =
- type$ "inbook" =
- or
- 'author.editor.full
- { type$ "proceedings" =
- 'editor.full
- 'author.full
- if$
- }
- if$
-}
-
-FUNCTION {output.bibitem}
-{ newline$
- "\bibitem[" write$
- label write$
- ")" make.full.names duplicate$ short.list =
- { pop$ }
- { * }
- if$
- "]{" * write$
- cite$ write$
- "}" write$
- newline$
- ""
- before.all 'output.state :=
-}
-
FUNCTION {n.dashify}
{ 't :=
""
{ t empty$ not }
{ t #1 #1 substring$ "-" =
{ t #1 #2 substring$ "--" = not
{ "--" *
t #2 global.max$ substring$ 't :=
}
{ { t #1 #1 substring$ "-" = }
{ "-" *
t #2 global.max$ substring$ 't :=
}
while$
}
if$
}
{ t #1 #1 substring$ *
t #2 global.max$ substring$ 't :=
}
if$
}
while$
}
-FUNCTION {format.date}
-{ year duplicate$ empty$
- { "empty year in " cite$ * warning$
- pop$ "" }
- 'skip$
- if$
- month empty$
- 'skip$
- { month
- " " * swap$ *
- }
- if$
- extra.label *
-}
-
FUNCTION {format.btitle}
{ title emphasize
}
FUNCTION {tie.or.space.connect}
{ duplicate$ text.length$ #3 <
{ "~" }
{ " " }
if$
swap$ * *
}
FUNCTION {either.or.check}
{ empty$
'pop$
{ "can't use both " swap$ * " fields in " * cite$ * warning$ }
if$
}
FUNCTION {format.bvolume}
{ volume empty$
{ "" }
{ "volume" volume tie.or.space.connect
series empty$
'skip$
{ " of " * series emphasize * }
if$
"volume and number" number either.or.check
}
if$
}
FUNCTION {format.number.series}
{ volume empty$
{ number empty$
{ series field.or.null }
{ output.state mid.sentence =
{ "number" }
{ "Number" }
if$
number tie.or.space.connect
series empty$
{ "there's a number but no series in " cite$ * warning$ }
{ " in " * series * }
if$
}
if$
}
{ "" }
if$
}
FUNCTION {format.edition}
{ edition empty$
{ "" }
{ output.state mid.sentence =
{ edition "l" change.case$ " edition" * }
{ edition "t" change.case$ " edition" * }
if$
}
if$
}
INTEGERS { multiresult }
FUNCTION {multi.page.check}
{ 't :=
#0 'multiresult :=
{ multiresult not
t empty$ not
and
}
{ t #1 #1 substring$
duplicate$ "-" =
swap$ duplicate$ "," =
swap$ "+" =
or or
{ #1 'multiresult := }
{ t #2 global.max$ substring$ 't := }
if$
}
while$
multiresult
}
FUNCTION {format.pages}
{ pages empty$
{ "" }
{ pages multi.page.check
{ "pages" pages n.dashify tie.or.space.connect }
{ "page" pages tie.or.space.connect }
if$
}
if$
}
-FUNCTION {format.eid}
-{ eid empty$
- { "" }
- { "art." eid tie.or.space.connect }
- if$
-}
-
FUNCTION {format.vol.num.pages}
{ volume field.or.null
number empty$
'skip$
- { "\penalty0 (" number * ")" * *
+ { "(" number * ")" * *
volume empty$
{ "there's a number but no volume in " cite$ * warning$ }
'skip$
if$
}
if$
pages empty$
'skip$
{ duplicate$ empty$
{ pop$ format.pages }
- { ":\penalty0 " * pages n.dashify * }
- if$
- }
- if$
-}
-
-FUNCTION {format.vol.num.eid}
-{ volume field.or.null
- number empty$
- 'skip$
- { "\penalty0 (" number * ")" * *
- volume empty$
- { "there's a number but no volume in " cite$ * warning$ }
- 'skip$
- if$
- }
- if$
- eid empty$
- 'skip$
- { duplicate$ empty$
- { pop$ format.eid }
- { ":\penalty0 " * eid * }
+ { ":" * pages n.dashify * }
if$
}
if$
}
FUNCTION {format.chapter.pages}
{ chapter empty$
'format.pages
{ type empty$
{ "chapter" }
{ type "l" change.case$ }
if$
chapter tie.or.space.connect
pages empty$
'skip$
{ ", " * format.pages * }
if$
}
if$
}
FUNCTION {format.in.ed.booktitle}
{ booktitle empty$
{ "" }
{ editor empty$
{ "In " booktitle emphasize * }
{ "In " format.editors * ", " * booktitle emphasize * }
if$
}
if$
}
-FUNCTION {empty.misc.check}
-{ author empty$ title empty$ howpublished empty$
- month empty$ year empty$ note empty$
- and and and and and
- key empty$ not and
- { "all relevant fields are empty in " cite$ * warning$ }
- 'skip$
- if$
-}
-
FUNCTION {format.thesis.type}
{ type empty$
'skip$
{ pop$
type "t" change.case$
}
if$
}
FUNCTION {format.tr.number}
{ type empty$
{ "Technical Report" }
'type
if$
number empty$
{ "t" change.case$ }
{ number tie.or.space.connect }
if$
}
FUNCTION {format.article.crossref}
-{ key empty$
- { journal empty$
- { "need key or journal for " cite$ * " to crossref " * crossref *
- warning$
- ""
- }
- { "In \emph{" journal * "}" * }
- if$
- }
- { "In " }
- if$
- " \citet{" * crossref * "}" *
+{ "In" % this is for apalike
+ " \cite{" * crossref * "}" *
}
FUNCTION {format.book.crossref}
{ volume empty$
{ "empty volume in " cite$ * "'s crossref of " * crossref * warning$
"In "
}
{ "Volume" volume tie.or.space.connect
" of " *
}
if$
- editor empty$
- editor field.or.null author field.or.null =
- or
- { key empty$
- { series empty$
- { "need editor, key, or series for " cite$ * " to crossref " *
- crossref * warning$
- "" *
- }
- { "\emph{" * series * "}" * }
- if$
- }
- 'skip$
- if$
- }
- 'skip$
- if$
- " \citet{" * crossref * "}" *
+ "\cite{" * crossref * "}" * % this is for apalike
}
FUNCTION {format.incoll.inproc.crossref}
-{ editor empty$
- editor field.or.null author field.or.null =
- or
- { key empty$
- { booktitle empty$
- { "need editor, key, or booktitle for " cite$ * " to crossref " *
- crossref * warning$
- ""
- }
- { "In \emph{" booktitle * "}" * }
- if$
- }
- { "In " }
- if$
- }
- { "In " }
- if$
- " \citet{" * crossref * "}" *
+{ "In" % this is for apalike
+ " \cite{" * crossref * "}" *
}
FUNCTION {article}
{ output.bibitem
format.authors "author" output.check
- author format.key output
+ author format.key output % special for
+ output.year.check % apalike
new.block
format.title "title" output.check
new.block
crossref missing$
{ journal emphasize "journal" output.check
- eid empty$
- { format.vol.num.pages output }
- { format.vol.num.eid output }
- if$
- format.date "year" output.check
+ format.vol.num.pages output
}
{ format.article.crossref output.nonnull
- eid empty$
- { format.pages output }
- { format.eid output }
- if$
+ format.pages output
}
if$
- format.issn output
- format.doi output
- format.eprint output
- format.url output
new.block
note output
fin.entry
}
FUNCTION {book}
{ output.bibitem
author empty$
{ format.editors "author and editor" output.check
editor format.key output
}
{ format.authors output.nonnull
crossref missing$
{ "author and editor" editor either.or.check }
'skip$
if$
}
if$
+ output.year.check % special for apalike
new.block
format.btitle "title" output.check
crossref missing$
{ format.bvolume output
new.block
format.number.series output
new.sentence
publisher "publisher" output.check
address output
}
{ new.block
format.book.crossref output.nonnull
}
if$
format.edition output
- format.date "year" output.check
- format.isbn output
- format.doi output
- format.eprint output
- format.url output
new.block
note output
fin.entry
}
FUNCTION {booklet}
{ output.bibitem
format.authors output
- author format.key output
+ author format.key output % special for
+ output.year.check % apalike
new.block
format.title "title" output.check
- howpublished address new.block.checkb
+ new.block
howpublished output
address output
- format.date output
- format.isbn output
- format.doi output
- format.eprint output
- format.url output
new.block
note output
fin.entry
}
FUNCTION {inbook}
{ output.bibitem
author empty$
{ format.editors "author and editor" output.check
editor format.key output
}
{ format.authors output.nonnull
crossref missing$
{ "author and editor" editor either.or.check }
'skip$
if$
}
if$
+ output.year.check % special for apalike
new.block
format.btitle "title" output.check
crossref missing$
{ format.bvolume output
format.chapter.pages "chapter and pages" output.check
new.block
format.number.series output
new.sentence
publisher "publisher" output.check
address output
}
{ format.chapter.pages "chapter and pages" output.check
new.block
format.book.crossref output.nonnull
}
if$
format.edition output
- format.date "year" output.check
- format.isbn output
- format.doi output
- format.eprint output
- format.url output
new.block
note output
fin.entry
}
FUNCTION {incollection}
{ output.bibitem
format.authors "author" output.check
- author format.key output
+ author format.key output % special for
+ output.year.check % apalike
new.block
format.title "title" output.check
new.block
crossref missing$
{ format.in.ed.booktitle "booktitle" output.check
format.bvolume output
format.number.series output
format.chapter.pages output
new.sentence
publisher "publisher" output.check
address output
format.edition output
- format.date "year" output.check
}
{ format.incoll.inproc.crossref output.nonnull
format.chapter.pages output
}
if$
- format.isbn output
- format.doi output
- format.eprint output
- format.url output
new.block
note output
fin.entry
}
FUNCTION {inproceedings}
{ output.bibitem
format.authors "author" output.check
- author format.key output
+ author format.key output % special for
+ output.year.check % apalike
new.block
format.title "title" output.check
new.block
crossref missing$
{ format.in.ed.booktitle "booktitle" output.check
format.bvolume output
format.number.series output
format.pages output
- address empty$
- { organization publisher new.sentence.checkb
- organization output
- publisher output
- format.date "year" output.check
- }
- { address output.nonnull
- format.date "year" output.check
- new.sentence
- organization output
- publisher output
- }
- if$
+ address output % for apalike
+ new.sentence % there's no year
+ organization output % here so things
+ publisher output % are simpler
}
{ format.incoll.inproc.crossref output.nonnull
format.pages output
}
if$
- format.isbn output
- format.doi output
- format.eprint output
- format.url output
new.block
note output
fin.entry
}
FUNCTION {conference} { inproceedings }
FUNCTION {manual}
{ output.bibitem
format.authors output
- author format.key output
+ author format.key output % special for
+ output.year.check % apalike
new.block
format.btitle "title" output.check
organization address new.block.checkb
organization output
address output
format.edition output
- format.date output
- format.url output
new.block
note output
fin.entry
}
FUNCTION {mastersthesis}
{ output.bibitem
format.authors "author" output.check
- author format.key output
+ author format.key output % special for
+ output.year.check % apalike
new.block
format.title "title" output.check
new.block
"Master's thesis" format.thesis.type output.nonnull
school "school" output.check
address output
- format.date "year" output.check
- format.url output
new.block
note output
fin.entry
}
FUNCTION {misc}
{ output.bibitem
format.authors output
- author format.key output
- title howpublished new.block.checkb
+ author format.key output % special for
+ output.year.check % apalike
+ new.block
format.title output
- howpublished new.block.checka
+ new.block
howpublished output
- format.date output
- format.issn output
- format.url output
new.block
note output
fin.entry
- empty.misc.check
}
FUNCTION {phdthesis}
{ output.bibitem
format.authors "author" output.check
- author format.key output
+ author format.key output % special for
+ output.year.check % apalike
new.block
format.btitle "title" output.check
new.block
"PhD thesis" format.thesis.type output.nonnull
school "school" output.check
address output
- format.date "year" output.check
- format.url output
new.block
note output
fin.entry
}
FUNCTION {proceedings}
{ output.bibitem
format.editors output
- editor format.key output
+ editor format.key output % special for
+ output.year.check % apalike
new.block
format.btitle "title" output.check
format.bvolume output
format.number.series output
- address output
- format.date "year" output.check
- new.sentence
- organization output
- publisher output
- format.isbn output
- format.doi output
- format.eprint output
- format.url output
+ address output % for apalike
+ new.sentence % we always output
+ organization output % a nonempty organization
+ publisher output % here
new.block
note output
fin.entry
}
FUNCTION {techreport}
{ output.bibitem
format.authors "author" output.check
- author format.key output
+ author format.key output % special for
+ output.year.check % apalike
new.block
format.title "title" output.check
new.block
format.tr.number output.nonnull
institution "institution" output.check
address output
- format.date "year" output.check
- format.url output
new.block
note output
fin.entry
}
FUNCTION {unpublished}
{ output.bibitem
format.authors "author" output.check
- author format.key output
+ author format.key output % special for
+ output.year.check % apalike
new.block
format.title "title" output.check
new.block
note "note" output.check
- format.date output
- format.url output
fin.entry
}
FUNCTION {default.type} { misc }
-
MACRO {jan} {"January"}
MACRO {feb} {"February"}
MACRO {mar} {"March"}
MACRO {apr} {"April"}
MACRO {may} {"May"}
MACRO {jun} {"June"}
MACRO {jul} {"July"}
MACRO {aug} {"August"}
MACRO {sep} {"September"}
MACRO {oct} {"October"}
MACRO {nov} {"November"}
MACRO {dec} {"December"}
-
-
MACRO {acmcs} {"ACM Computing Surveys"}
MACRO {acta} {"Acta Informatica"}
MACRO {cacm} {"Communications of the ACM"}
MACRO {ibmjrd} {"IBM Journal of Research and Development"}
MACRO {ibmsj} {"IBM Systems Journal"}
MACRO {ieeese} {"IEEE Transactions on Software Engineering"}
MACRO {ieeetc} {"IEEE Transactions on Computers"}
MACRO {ieeetcad}
{"IEEE Transactions on Computer-Aided Design of Integrated Circuits"}
MACRO {ipl} {"Information Processing Letters"}
MACRO {jacm} {"Journal of the ACM"}
MACRO {jcss} {"Journal of Computer and System Sciences"}
MACRO {scp} {"Science of Computer Programming"}
MACRO {sicomp} {"SIAM Journal on Computing"}
MACRO {tocs} {"ACM Transactions on Computer Systems"}
MACRO {tods} {"ACM Transactions on Database Systems"}
MACRO {tog} {"ACM Transactions on Graphics"}
MACRO {toms} {"ACM Transactions on Mathematical Software"}
MACRO {toois} {"ACM Transactions on Office Information Systems"}
MACRO {toplas} {"ACM Transactions on Programming Languages and Systems"}
MACRO {tcs} {"Theoretical Computer Science"}
-
READ
FUNCTION {sortify}
{ purify$
"l" change.case$
}
INTEGERS { len }
FUNCTION {chop.word}
{ 's :=
'len :=
s #1 len substring$ =
{ s len #1 + global.max$ substring$ }
's
if$
}
+% There are three apalike cases: one person (Jones),
+% two (Jones and de~Bruijn), and more (Jones et~al.).
+% This function is much like format.crossref.editors.
+%
FUNCTION {format.lab.names}
{ 's :=
s #1 "{vv~}{ll}" format.name$
s num.names$ duplicate$
#2 >
{ pop$ " et~al." * }
{ #2 <
'skip$
{ s #2 "{ff }{vv }{ll}{ jj}" format.name$ "others" =
{ " et~al." * }
{ " and " * s #2 "{vv~}{ll}" format.name$ * }
if$
}
if$
}
if$
}
FUNCTION {author.key.label}
{ author empty$
{ key empty$
{ cite$ #1 #3 substring$ }
- 'key
+ 'key % apalike uses the whole key
if$
}
{ author format.lab.names }
if$
}
FUNCTION {author.editor.key.label}
{ author empty$
{ editor empty$
{ key empty$
{ cite$ #1 #3 substring$ }
- 'key
+ 'key % apalike uses the whole key
if$
}
{ editor format.lab.names }
if$
}
{ author format.lab.names }
if$
}
-FUNCTION {author.key.organization.label}
-{ author empty$
- { key empty$
- { organization empty$
- { cite$ #1 #3 substring$ }
- { "The " #4 organization chop.word #3 text.prefix$ }
- if$
- }
- 'key
- if$
- }
- { author format.lab.names }
- if$
-}
-
-FUNCTION {editor.key.organization.label}
+FUNCTION {editor.key.label}
{ editor empty$
{ key empty$
- { organization empty$
- { cite$ #1 #3 substring$ }
- { "The " #4 organization chop.word #3 text.prefix$ }
- if$
- }
- 'key
+ { cite$ #1 #3 substring$ }
+ 'key % apalike uses the whole key, no organization
if$
}
{ editor format.lab.names }
if$
}
-FUNCTION {calc.short.authors}
+FUNCTION {calc.label}
{ type$ "book" =
type$ "inbook" =
or
'author.editor.key.label
{ type$ "proceedings" =
- 'editor.key.organization.label
- { type$ "manual" =
- 'author.key.organization.label
- 'author.key.label
- if$
- }
+ 'editor.key.label % apalike ignores organization
+ 'author.key.label % for labeling and sorting
if$
}
if$
- 'short.list :=
-}
-
-FUNCTION {calc.label}
-{ calc.short.authors
- short.list
- "("
- *
- year duplicate$ empty$
- short.list key field.or.null = or
- { pop$ "" }
- 'skip$
- if$
+ ", " % these three lines are
+ * % for apalike, which
+ year field.or.null purify$ #-1 #4 substring$ % uses all four digits
*
'label :=
}
-INTEGERS { seq.num }
+FUNCTION {sort.format.names}
+{ 's :=
+ #1 'nameptr :=
+ ""
+ s num.names$ 'numnames :=
+ numnames 'namesleft :=
+ { namesleft #0 > }
+ { nameptr #1 >
+ { " " * }
+ 'skip$
+ if$ % apalike uses initials
+ s nameptr "{vv{ } }{ll{ }}{ f{ }}{ jj{ }}" format.name$ 't := % <= here
+ nameptr numnames = t "others" = and
+ { "et al" * }
+ { t sortify * }
+ if$
+ nameptr #1 + 'nameptr :=
+ namesleft #1 - 'namesleft :=
+ }
+ while$
+}
-FUNCTION {init.seq}
-{ #0 'seq.num :=}
+FUNCTION {sort.format.title}
+{ 't :=
+ "A " #2
+ "An " #3
+ "The " #4 t chop.word
+ chop.word
+ chop.word
+ sortify
+ #1 global.max$ substring$
+}
-EXECUTE {init.seq}
+FUNCTION {author.sort}
+{ author empty$
+ { key empty$
+ { "to sort, need author or key in " cite$ * warning$
+ ""
+ }
+ { key sortify }
+ if$
+ }
+ { author sort.format.names }
+ if$
+}
-FUNCTION {int.to.fix}
-{ "000000000" swap$ int.to.str$ *
- #-1 #10 substring$
+FUNCTION {author.editor.sort}
+{ author empty$
+ { editor empty$
+ { key empty$
+ { "to sort, need author, editor, or key in " cite$ * warning$
+ ""
+ }
+ { key sortify }
+ if$
+ }
+ { editor sort.format.names }
+ if$
+ }
+ { author sort.format.names }
+ if$
}
+FUNCTION {editor.sort}
+{ editor empty$
+ { key empty$
+ { "to sort, need editor or key in " cite$ * warning$
+ ""
+ }
+ { key sortify }
+ if$
+ }
+ { editor sort.format.names }
+ if$
+}
+% apalike uses two sorting passes; the first one sets the
+% labels so that the `a's, `b's, etc. can be computed;
+% the second pass puts the references in "correct" order.
+% The presort function is for the first pass. It computes
+% label, sort.label, and title, and then concatenates.
FUNCTION {presort}
{ calc.label
label sortify
" "
*
- seq.num #1 + 'seq.num :=
- seq.num int.to.fix
- 'sort.label :=
- sort.label *
+ type$ "book" =
+ type$ "inbook" =
+ or
+ 'author.editor.sort
+ { type$ "proceedings" =
+ 'editor.sort
+ 'author.sort
+ if$
+ }
+ if$
+ #1 entry.max$ substring$ % for
+ 'sort.label := % apalike
+ sort.label % style
+ *
+ " "
+ *
+ title field.or.null
+ sort.format.title
+ *
#1 entry.max$ substring$
'sort.key$ :=
}
ITERATE {presort}
-SORT
+SORT % by label, sort.label, title---for final label calculation
-STRINGS { longest.label last.label next.extra }
+STRINGS { last.label next.extra } % apalike labels are only for the text;
-INTEGERS { longest.label.width last.extra.num number.label }
+INTEGERS { last.extra.num } % there are none in the bibliography
-FUNCTION {initialize.longest.label}
-{ "" 'longest.label :=
- #0 int.to.chr$ 'last.label :=
+FUNCTION {initialize.extra.label.stuff} % and hence there is no `longest.label'
+{ #0 int.to.chr$ 'last.label :=
"" 'next.extra :=
- #0 'longest.label.width :=
#0 'last.extra.num :=
- #0 'number.label :=
}
FUNCTION {forward.pass}
{ last.label label =
{ last.extra.num #1 + 'last.extra.num :=
last.extra.num int.to.chr$ 'extra.label :=
}
{ "a" chr.to.int$ 'last.extra.num :=
"" 'extra.label :=
label 'last.label :=
}
if$
- number.label #1 + 'number.label :=
}
FUNCTION {reverse.pass}
{ next.extra "b" =
{ "a" 'extra.label := }
'skip$
if$
- extra.label 'next.extra :=
- extra.label
- duplicate$ empty$
- 'skip$
- { "{\natexlab{" swap$ * "}}" * }
- if$
- 'extra.label :=
label extra.label * 'label :=
+ extra.label 'next.extra :=
}
-EXECUTE {initialize.longest.label}
+EXECUTE {initialize.extra.label.stuff}
ITERATE {forward.pass}
REVERSE {reverse.pass}
+% Now that the label is right we sort for real,
+% on sort.label then year then title. This is
+% for the second sorting pass.
FUNCTION {bib.sort.order}
-{ sort.label 'sort.key$ :=
+{ sort.label
+ " "
+ *
+ year field.or.null sortify
+ *
+ " "
+ *
+ title field.or.null
+ sort.format.title
+ *
+ #1 entry.max$ substring$
+ 'sort.key$ :=
}
ITERATE {bib.sort.order}
-SORT
+SORT % by sort.label, year, title---giving final bibliography order
FUNCTION {begin.bib}
-{ preamble$ empty$
+{ preamble$ empty$ % no \etalchar in apalike
'skip$
{ preamble$ write$ newline$ }
if$
- "\begin{thebibliography}{" number.label int.to.str$ * "}" *
- write$ newline$
- "\providecommand{\natexlab}[1]{#1}"
- write$ newline$
- "\providecommand{\url}[1]{\texttt{#1}}"
- write$ newline$
- "\expandafter\ifx\csname urlstyle\endcsname\relax"
- write$ newline$
- " \providecommand{\eprint}[1]{eprint: #1}\else"
- write$ newline$
- " \providecommand{\eprint}{eprint: \begingroup \urlstyle{rm}\Url}\fi"
- write$ newline$
- "\expandafter\ifx\csname urlstyle\endcsname\relax"
- write$ newline$
- " \providecommand{\doi}[1]{doi: #1}\else"
- write$ newline$
- " \providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi"
- write$ newline$
+ "\begin{thebibliography}{}" write$ newline$ % no labels in apalike
}
EXECUTE {begin.bib}
EXECUTE {init.state.consts}
ITERATE {call.type$}
FUNCTION {end.bib}
{ newline$
"\end{thebibliography}" write$ newline$
}
EXECUTE {end.bib}
diff --git a/docs/theory/examples/axisymmetric_linear-elasticity.py b/docs/theory/examples/axisymmetric_linear-elasticity.py
index 0af7a2e..ad4a9c8 100644
--- a/docs/theory/examples/axisymmetric_linear-elasticity.py
+++ b/docs/theory/examples/axisymmetric_linear-elasticity.py
@@ -1,310 +1,406 @@
+# ==================================================================================================
+#
+# (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
+#
+# ==================================================================================================
import numpy as np
np.set_printoptions(linewidth=200)
+# ==================================================================================================
+# tensor products
+# ==================================================================================================
+
+ddot22 = lambda A2,B2: np.einsum('...ij ,...ji->... ',A2,B2)
+ddot42 = lambda A4,B2: np.einsum('...ijkl,...lk->...ij ',A4,B2)
+dyad22 = lambda A2,B2: np.einsum('...ij ,...kl->...ijkl',A2,B2)
+ddot43 = lambda A4,B3: np.einsum('...ijkl,...lkm->...ijm',A4,B3)
+ddot33 = lambda A3,B3: np.einsum('...ijk ,...kjl->...il ',A3,B3)
+dot31 = lambda A3,B1: np.einsum('...ijk ,...k ->...ij ',A3,B1)
+transpose3 = lambda A3 : np.einsum('...ijk ->...kji',A3 )
+
# ==================================================================================================
# mesh definition
# ==================================================================================================
# number of elements
nz = 20
nr = 16
# mesh dimensions
nelem = nz * nr # number of elements
nnode = (nz+1) * (nr+1) # number of nodes
nne = 4 # number of nodes per element
-nd = 2 # number of dimensions
-ndof = nnode * nd # total number of degrees of freedom
+ndim = 2 # number of dimensions
+ndof = nnode * ndim # total number of degrees of freedom
-# coordinates and connectivity: zero-initialize
-coor = np.zeros((nnode,nd ), dtype='float')
-conn = np.zeros((nelem,nne), dtype='int' )
+# zero-initialise coordinates, displacements, and connectivity
+coor = np.zeros((nnode,ndim), dtype='float')
+disp = np.zeros((nnode,ndim), dtype='float')
+conn = np.zeros((nelem,nne ), dtype='int' )
-# coordinates: set
+# coordinates
# - grid of points
z,r = np.meshgrid(np.linspace(0,1,nz+1), np.linspace(0,1,nr+1))
# - store from grid of points
coor[:,0] = z.ravel()
coor[:,1] = r.ravel()
-# connectivity: set
+# connectivity
# - grid of node numbers
inode = np.arange(nnode).reshape(nr+1, nz+1)
# - store from grid of node numbers
conn[:,0] = inode[ :-1, :-1].ravel()
conn[:,1] = inode[ :-1,1: ].ravel()
conn[:,2] = inode[1: ,1: ].ravel()
conn[:,3] = inode[1: , :-1].ravel()
+# - node sets
+nodesLeft = inode[ :, 0]
+nodesRight = inode[ :,-1]
+nodesBottom = inode[ 0, :]
+nodesTop = inode[-1, :]
# DOF-numbers per node
-dofs = np.arange(nnode*nd).reshape(nnode,nd)
+dofs = np.arange(nnode*ndim).reshape(nnode,ndim)
# ==================================================================================================
# quadrature definition
# ==================================================================================================
# integration point coordinates (local element coordinates)
Xi = np.array([
[-1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), +1./np.sqrt(3.)],
[-1./np.sqrt(3.), +1./np.sqrt(3.)],
])
# integration point weights
W = np.array([
[1.],
[1.],
[1.],
[1.],
])
# number of integration points
nip = 4
# ==================================================================================================
-# integration point tensors and operations
+# B-matrix at each integration point
# ==================================================================================================
-# tensors products
-ddot22 = lambda A2,B2: np.einsum('...ij ,...ji->... ',A2,B2)
-ddot42 = lambda A4,B2: np.einsum('...ijkl,...lk->...ij ',A4,B2)
-dyad22 = lambda A2,B2: np.einsum('...ij ,...kl->...ijkl',A2,B2)
-ddot43 = lambda A4,B3: np.einsum('...ijkl,...lkm->...ijm',A4,B3)
-ddot33 = lambda A3,B3: np.einsum('...ijk ,...kjl->...il ',A3,B3)
-dot31 = lambda A3,B1: np.einsum('...ijk ,...k ->...ij ',A3,B1)
-transpose3 = lambda A3 : np.einsum('...ijk ->...kji',A3 )
-
-# identity tensors (per integration point)
-i = np.eye(3)
-I = np.einsum('ij ,...' , i ,np.ones([nelem,nip]))
-I4 = np.einsum('ijkl,...->...ijkl',np.einsum('il,jk',i,i),np.ones([nelem,nip]))
-I4rt = np.einsum('ijkl,...->...ijkl',np.einsum('ik,jl',i,i),np.ones([nelem,nip]))
-I4s = (I4+I4rt)/2.
-II = dyad22(I,I)
-I4d = I4s-II/3.
-
-# bulk and shear modulus
-K = 0.8333333333333333
-G = 0.3846153846153846
-
-# elasticity tensor (per integration point)
-C4 = K * II + 2. * G * I4d
-
-# ==================================================================================================
-# B-matrices at each integration point
-# ==================================================================================================
-
-# allocate matrix
-B = np.zeros((nelem,nip,nne,3,3,3))
-vol = np.zeros((nelem,nip))
+# allocate
+B = np.empty((nelem,nip,nne,3,3,3))
+vol = np.empty((nelem,nip))
# loop over nodes
for e, nodes in enumerate(conn):
# nodal coordinates
xe = coor[nodes,:]
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# shape functions
N = np.array([
.25 * (1.-xi[0]) * (1.-xi[1]),
.25 * (1.+xi[0]) * (1.-xi[1]),
.25 * (1.+xi[0]) * (1.+xi[1]),
.25 * (1.-xi[0]) * (1.+xi[1]),
])
# shape function gradients (w.r.t. the local element coordinates)
dNdxi = np.array([
[-.25*(1.-xi[1]), -.25*(1.-xi[0])],
[+.25*(1.-xi[1]), -.25*(1.+xi[0])],
[+.25*(1.+xi[1]), +.25*(1.+xi[0])],
[-.25*(1.+xi[1]), +.25*(1.-xi[0])],
])
# Jacobian
- J = np.einsum('mi,mj->ij', dNdxi, xe)
- Jinv = np.linalg.inv(J)
- Jdet = np.linalg.det(J)
+ Je = np.einsum('mi,mj->ij', dNdxi, xe)
+ invJe = np.linalg.inv(Je)
+ detJe = np.linalg.det(Je)
# shape function gradients (w.r.t. the global coordinates)
- dNdx = np.einsum('ij,mj->mi', Jinv, dNdxi)
+ dNdxe = np.einsum('ij,mj->mi', invJe, dNdxi)
- # global coordinates of the integration point
+ # global coordinates of the integration point, extract the radius
xq = np.einsum('m,mi->i', N, xe)
rq = xq[1]
# compute B-matrix and integration-point volume and store for later use
for m in range(nne):
- Bm = np.zeros((3,3,3))
+ Be = np.zeros((3,3,3))
- Bm[0,1,1] = -1./rq * N[m] # B(m, r \theta \theta )
- Bm[1,1,0] = +1./rq * N[m] # B(m, \theta \theta r )
+ Be[0,0,0] = dNdxe[m,1] # B(m, r r r )
+ Be[0,2,2] = dNdxe[m,1] # B(m, r z z )
+ Be[1,1,0] = +1./rq * N[m] # B(m, \theta \theta r )
+ Be[2,0,0] = dNdxe[m,0] # B(m, z r r )
+ Be[2,2,2] = dNdxe[m,0] # B(m, z z z )
- Bm[0,0,0] = dNdx[m,1] # B(m, r r r )
- Bm[0,1,1] = dNdx[m,1] # B(m, r \theta \theta )
- Bm[0,2,2] = dNdx[m,1] # B(m, r z z )
+ B[e,q,m,:,:,:] = Be
- Bm[2,0,0] = dNdx[m,0] # B(m, z r r )
- Bm[2,1,1] = dNdx[m,0] # B(m, z \theta \theta )
- Bm[2,2,2] = dNdx[m,0] # B(m, z z z )
+ vol[e,q] = w * detJe * rq * 2. * np.pi
- B[e,q,m,:,:,:] = Bm
+# ==================================================================================================
+# stiffness tensor at each integration point (provides constitutive response and 'tangent')
+# ==================================================================================================
- vol[e,q] = w * Jdet * rq * 2. * np.pi
+# identity tensors (per integration point)
+i = np.eye(3)
+I = np.einsum('ij ,...' , i ,np.ones([nelem,nip]))
+I4 = np.einsum('ijkl,...->...ijkl',np.einsum('il,jk',i,i),np.ones([nelem,nip]))
+I4rt = np.einsum('ijkl,...->...ijkl',np.einsum('ik,jl',i,i),np.ones([nelem,nip]))
+I4s = (I4+I4rt)/2.
+II = dyad22(I,I)
+I4d = I4s-II/3.
+
+# bulk and shear modulus (homogeneous)
+K = 0.8333333333333333 * np.ones((nelem, nip))
+G = 0.3846153846153846 * np.ones((nelem, nip))
+
+# elasticity tensor (per integration point)
+C4 = np.einsum('eq,eqijkl->eqijkl', K, II) + 2. * np.einsum('eq,eqijkl->eqijkl', G, I4d)
# ==================================================================================================
-# assemble stiffness matrix
+# strain from nodal displacements, stress from constitutive response
# ==================================================================================================
-# allocate matrix
-K = np.zeros((ndof, ndof))
+# allocate strain tensor per integration point
+Eps = np.empty((nelem,nip,3,3))
# loop over nodes
for e, nodes in enumerate(conn):
- # allocate element stiffness matrix
- Ke = np.zeros((nne*nd, nne*nd))
+ # nodal displacements in 3-d
+ # u_theta = 0
+ # (z,r) -> (r,theta,z)
+ ue = np.zeros((nne,3))
+ ue[:,0] = disp[nodes,1]
+ ue[:,2] = disp[nodes,0]
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
# alias integration point values
- C4q = C4[e,q,:,:,:,:]
- dV = vol[e,q]
+ Be = B[e,q,:,:,:,:]
- # assemble to element stiffness matrix
- for m in range(nne):
+ # displacement gradient
+ gradu = np.einsum('mijk,mk->ij', Be, ue)
+
+ # compute strain tensor, and store per integration point
+ Eps[e,q] = .5 * ( gradu + gradu.T )
+
+# constitutive response: strain tensor -> stress tensor (per integration point)
+Sig = ddot42(C4, Eps)
- Bm = B[e,q,m,:,:,:]
+# ==================================================================================================
+# internal force from stress
+# ==================================================================================================
- for n in range(nne):
+# allocate internal force
+fint = np.zeros((ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
- Bn = B[e,q,n,:,:,:]
+ # allocate element internal force
+ fe = np.zeros((nne*ndim))
- Kmn = ddot33(transpose3(Bm),ddot43(C4q,Bn))
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
- iim = np.array([ m*nd+0, m*nd+1 ])
- iin = np.array([ n*nd+0, n*nd+1 ])
+ # alias integration point values
+ sig = Sig[e,q,:,:]
+ Be = B [e,q,:,:,:,:]
+ dV = vol[e,q]
- Ke[np.ix_(iim,iin)] += Kmn[np.ix_([2,0], [2,0])] * dV
+ # assemble to element internal force
+ # Bm = B[e,q,m,:,:,:]
+ # fm = ddot32(transpose3(Bm), sig) * dV
+ # fe[[m*ndim+0, m*ndim+1]] = [fm[m,2], fm[m,0]]
+ fq = np.einsum('mijk,ij->mk', Be, sig) * dV
+ fe += fq[:,[2,0]].reshape(nne*ndim)
# assemble to global stiffness matrix
iie = dofs[nodes,:].ravel()
+ fint[iie] += fe
+
+# ==================================================================================================
+# stiffness matrix from 'tangent'
+# ==================================================================================================
+# allocate stiffness matrix
+K = np.zeros((ndof, ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element stiffness matrix
+ Ke = np.zeros((nne*ndim, nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ c4 = C4 [e,q,:,:,:,:]
+ Be = B [e,q,:,:,:,:]
+ dV = vol[e,q]
+
+ # assemble to element stiffness matrix
+ # Bm = B[e,q,m,:,:,:]
+ # Bn = B[e,q,n,:,:,:]
+ # Kmn = ddot33(transpose3(Bm),ddot43(C4,Bn)) * dV
+ # Ke[[m*ndim+0,m*ndim+1], [n*ndim+0,n*ndim+1]] += Kmn[[2,0], [2,0]]
+ Kq = np.einsum('mabc,abde,nedf->mcnf', Be, c4, Be) * dV
+ Ke += Kq[np.ix_(np.arange(nne),[2,0],np.arange(nne),[2,0])].reshape(nne*ndim, nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
K[np.ix_(iie,iie)] += Ke
# ==================================================================================================
# partition and solve
# ==================================================================================================
-# zero-initialize displacements and forces
-f = np.zeros((ndof))
-u = np.zeros((ndof))
+# prescribed external force: zero on all free DOFs
+# (other DOFs ignored in solution, the reaction forces on the DOFs are computed below)
+fext = np.zeros((ndof))
+
+# DOF-partitioning: ['u'nknown, 'p'rescribed]
+# - prescribed
+iip = np.hstack((
+ dofs[nodesBottom,1],
+ dofs[nodesLeft ,0],
+ dofs[nodesRight ,0],
+))
+# - unknown
+iiu = np.setdiff1d(dofs.ravel(), iip)
# fixed displacements
-# - zero-initialize
-iip = np.empty((0),dtype='int' )
-up = np.empty((0),dtype='float')
-# - 'r = 0' : 'u_r = 0'
-idof = dofs[inode[0,:], 1]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.0 * np.ones(idof.shape) ))
-# - 'z = 0' : 'u_z = 0'
-idof = dofs[inode[:,0], 0]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.0 * np.ones(idof.shape) ))
-# - 'z = 1' : 'u_z = 0.1'
-idof = dofs[inode[:,-1], 0]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.1 * np.ones(idof.shape) ))
-
-# free displacements
-iiu = np.setdiff1d(dofs.ravel(), iip)
+up = np.hstack((
+ 0.0 * np.ones((len(nodesBottom))),
+ 0.0 * np.ones((len(nodesLeft ))),
+ 0.1 * np.ones((len(nodesRight ))),
+))
+
+# residual force
+r = fext - fint
# partition
# - stiffness matrix
-Kuu = K[np.ix_(iiu, iiu)]
-Kup = K[np.ix_(iiu, iip)]
-Kpu = K[np.ix_(iip, iiu)]
-Kpp = K[np.ix_(iip, iip)]
-# - external force
-fu = f[iiu]
+Kuu = K[np.ix_(iiu, iiu)]
+Kup = K[np.ix_(iiu, iip)]
+Kpu = K[np.ix_(iip, iiu)]
+Kpp = K[np.ix_(iip, iip)]
+# - residual force
+ru = r[iiu]
-# solve
-uu = np.linalg.solve(Kuu, fu - Kup.dot(up))
-fp = Kpu.dot(uu) + Kpp.dot(up)
+# solve for unknown displacement DOFs
+uu = np.linalg.solve(Kuu, ru - Kup.dot(up))
# assemble
+u = np.empty((ndof))
u[iiu] = uu
u[iip] = up
-f[iip] = fp
-# convert to nodal displacement and forces
+# convert to nodal displacements
disp = u[dofs]
-fext = f[dofs]
# ==================================================================================================
-# compute strain and stress
+# strain from nodal displacements, stress from constitutive response
# ==================================================================================================
-# zero-initialize
-# - strain tensor per integration point
-eps = np.zeros((nelem,nip,3,3))
-# - nodal displacement in 3-d
-um = np.zeros((3))
+# allocate strain tensor per integration point
+Eps = np.empty((nelem,nip,3,3))
# loop over nodes
for e, nodes in enumerate(conn):
- # nodal displacements
- ue = disp[nodes,:]
+ # nodal displacements in 3-d
+ # u_theta = 0
+ # (z,r) -> (r,theta,z)
+ ue = np.zeros((nne,3))
+ ue[:,0] = disp[nodes,1]
+ ue[:,2] = disp[nodes,0]
# loop over integration points
for q, (w, xi) in enumerate(zip(W, Xi)):
- # zero-initialize displacements gradient
- gradu = np.zeros((3,3))
+ # alias integration point values
+ Be = B[e,q,:,:,:,:]
- # compute displacement gradient
- for m in range(nne):
+ # displacement gradient
+ gradu = np.einsum('mijk,mk->ij', Be, ue)
- # alias
- # - B-matrix
- Bm = B[e,q,m,:,:,:]
- # - nodal displacement in 3-d (theta-component always zero)
- um[np.ix_([2,0])] = ue[m,:]
+ # compute strain tensor, and store per integration point
+ Eps[e,q] = .5 * ( gradu + gradu.T )
- # update
- gradu += dot31(Bm, um)
+# constitutive response: strain tensor -> stress tensor (per integration point)
+Sig = ddot42(C4, Eps)
- # compute strain tensor, and store per integration point
- eps[e,q] = .5 * ( gradu + gradu.T )
+# ==================================================================================================
+# internal force from stress, reaction (external) force from internal force on fixed nodes
+# ==================================================================================================
-# compute stress tensor (per integration point)
-sig = ddot42(C4,eps)
+# allocate internal force
+fint = np.zeros((ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element internal force
+ fe = np.zeros((nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ sig = Sig[e,q,:,:]
+ Be = B [e,q,:,:,:,:]
+ dV = vol[e,q]
+
+ # assemble to element internal force
+ # Bm = B[e,q,m,:,:,:]
+ # fm = ddot32(transpose3(Bm), sig) * dV
+ # fe[[m*ndim+0, m*ndim+1]] = [fm[m,2], fm[m,0]]
+ fq = np.einsum('mijk,ij->mk', Be, sig) * dV
+ fe += fq[:,[2,0]].reshape(nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
+ fint[iie] += fe
+
+# reaction force
+fext[iip] = fint[iip]
# ==================================================================================================
# plot
# ==================================================================================================
import matplotlib.pyplot as plt
-fig, axes = plt.subplots(ncols=2, figsize=(16,8))
+fig, axes = plt.subplots(ncols=2, figsize=(16,8))
+
+# reconstruct external force as nodal vectors
+fext = fext[dofs]
+# plot original nodal positions + displacement field as arrows
axes[0].scatter(coor[:,0], coor[:,1])
axes[0].quiver (coor[:,0], coor[:,1], disp[:,0], disp[:,1])
+# plot new nodal positions + external (reaction) force field as arrows
axes[1].scatter((coor+disp)[:,0], (coor+disp)[:,1])
axes[1].quiver ((coor+disp)[:,0], (coor+disp)[:,1], fext[:,0], fext[:,1], scale=.4)
+# fix axes limits
for axis in axes:
axis.set_xlim([-0.4, 1.5])
axis.set_ylim([-0.4, 1.5])
plt.show()
diff --git a/docs/theory/examples/cartesian_linear-elasticity.py b/docs/theory/examples/cartesian_linear-elasticity.py
index 7ab7124..f0e11e1 100644
--- a/docs/theory/examples/cartesian_linear-elasticity.py
+++ b/docs/theory/examples/cartesian_linear-elasticity.py
@@ -1,239 +1,370 @@
-
-import numpy as np
-
-np.set_printoptions(linewidth=200)
-
# ==================================================================================================
-# identity tensors
+#
+# (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
+#
# ==================================================================================================
-II = np.zeros((3,3,3,3))
-I4 = np.zeros((3,3,3,3))
-I4rt = np.zeros((3,3,3,3))
-Is = np.zeros((3,3,3,3))
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == j and k == l ):
- II[i,j,k,l] = 1.
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == l and j == k ):
- I4[i,j,k,l] = 1.
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == k and j == l ):
- I4rt[i,j,k,l] = 1.
+import numpy as np
-I4s = ( I4 + I4rt ) / 2.
-I4d = ( I4s - II/3. )
+np.set_printoptions(linewidth=200)
# ==================================================================================================
-# elasticity tensor
+# tensor products
# ==================================================================================================
-# bulk and shear modulus
-K = 0.8333333333333333
-G = 0.3846153846153846
-
-# elasticity tensor (3d)
-C4 = K * II + 2. * G * I4d
+ddot22 = lambda A2,B2: np.einsum('...ij ,...ji->... ',A2,B2)
+ddot42 = lambda A4,B2: np.einsum('...ijkl,...lk->...ij ',A4,B2)
+dyad22 = lambda A2,B2: np.einsum('...ij ,...kl->...ijkl',A2,B2)
# ==================================================================================================
# mesh definition
# ==================================================================================================
# number of elements
-nx = 10
-ny = 10
+nx = 20
+ny = 20
# mesh dimensions
nelem = nx * ny # number of elements
nnode = (nx+1) * (ny+1) # number of nodes
nne = 4 # number of nodes per element
-nd = 2 # number of dimensions
-ndof = nnode * nd # total number of degrees of freedom
+ndim = 2 # number of dimensions
+ndof = nnode * ndim # total number of degrees of freedom
# out-of-plane thickness
thick = 1.
-# coordinates and connectivity: zero-initialize
-coor = np.zeros((nnode,nd ), dtype='float')
-conn = np.zeros((nelem,nne), dtype='int' )
+# zero-initialise coordinates, displacements, and connectivity
+coor = np.zeros((nnode,ndim), dtype='float')
+disp = np.zeros((nnode,ndim), dtype='float')
+conn = np.zeros((nelem,nne ), dtype='int' )
-# coordinates: set
+# coordinates
# - grid of points
x,y = np.meshgrid(np.linspace(0,1,nx+1), np.linspace(0,1,ny+1))
# - store from grid of points
coor[:,0] = x.ravel()
coor[:,1] = y.ravel()
-# connectivity: set
+# connectivity
# - grid of node numbers
inode = np.arange(nnode).reshape(ny+1, nx+1)
# - store from grid of node numbers
conn[:,0] = inode[ :-1, :-1].ravel()
conn[:,1] = inode[ :-1,1: ].ravel()
conn[:,2] = inode[1: ,1: ].ravel()
conn[:,3] = inode[1: , :-1].ravel()
+# - node sets
+nodesLeft = inode[ :, 0]
+nodesRight = inode[ :,-1]
+nodesBottom = inode[ 0, :]
+nodesTop = inode[-1, :]
# DOF-numbers per node
-dofs = np.arange(nnode*nd).reshape(nnode,nd)
+dofs = np.arange(nnode*ndim).reshape(nnode,ndim)
# ==================================================================================================
# quadrature definition
# ==================================================================================================
# integration point coordinates (local element coordinates)
Xi = np.array([
[-1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), +1./np.sqrt(3.)],
[-1./np.sqrt(3.), +1./np.sqrt(3.)],
])
# integration point weights
W = np.array([
[1.],
[1.],
[1.],
[1.],
])
+# number of integration points
+nip = 4
+
# ==================================================================================================
-# assemble stiffness matrix
+# gradient of the shape functions at each integration point
# ==================================================================================================
-# allocate matrix
-K = np.zeros((ndof, ndof))
+# allocate
+dNx = np.empty((nelem,nip,nne,ndim))
+vol = np.empty((nelem,nip))
# loop over nodes
-for e in conn:
+for e, nodes in enumerate(conn):
- # - nodal coordinates
- xe = coor[e,:]
+ # nodal coordinates
+ xe = coor[nodes,:]
- # - allocate element stiffness matrix
- Ke = np.zeros((nne*nd, nne*nd))
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
- # - numerical quadrature
- for w, xi in zip(W, Xi):
-
- # -- shape function gradients (w.r.t. the local element coordinates)
+ # shape function gradients (w.r.t. the local element coordinates)
dNdxi = np.array([
[-.25*(1.-xi[1]), -.25*(1.-xi[0])],
[+.25*(1.-xi[1]), -.25*(1.+xi[0])],
[+.25*(1.+xi[1]), +.25*(1.+xi[0])],
[-.25*(1.+xi[1]), +.25*(1.-xi[0])],
])
- # -- Jacobian
- J = np.zeros((nd, nd))
+ # Jacobian
+ Je = np.einsum('mi,mj->ij', dNdxi, xe)
+ invJe = np.linalg.inv(Je)
+ detJe = np.linalg.det(Je)
+
+ # shape function gradients (w.r.t. the global coordinates)
+ dNdxe = np.einsum('ij,mj->mi', invJe, dNdxi)
+
+ # store for later use
+ dNx[e,q,:,:] = dNdxe
+ vol[e,q] = w * detJe * thick
+
+# ==================================================================================================
+# stiffness tensor at each integration point (provides constitutive response and 'tangent')
+# ==================================================================================================
+
+# identity tensors (per integration point)
+i = np.eye(3)
+I = np.einsum('ij ,...' , i ,np.ones([nelem,nip]))
+I4 = np.einsum('ijkl,...->...ijkl',np.einsum('il,jk',i,i),np.ones([nelem,nip]))
+I4rt = np.einsum('ijkl,...->...ijkl',np.einsum('ik,jl',i,i),np.ones([nelem,nip]))
+I4s = (I4+I4rt)/2.
+II = dyad22(I,I)
+I4d = I4s-II/3.
- for m in range(nne):
- for i in range(nd):
- for j in range(nd):
- J[i,j] += dNdxi[m,i] * xe[m,j]
+# bulk and shear modulus (homogeneous)
+K = 0.8333333333333333 * np.ones((nelem, nip))
+G = 0.3846153846153846 * np.ones((nelem, nip))
+
+# elasticity tensor (per integration point)
+C4 = np.einsum('eq,eqijkl->eqijkl', K, II) + 2. * np.einsum('eq,eqijkl->eqijkl', G, I4d)
+
+# ==================================================================================================
+# strain from nodal displacements, stress from constitutive response
+# ==================================================================================================
+
+# zero-initialise strain tensor per integration point
+# (plain strain -> all z-components are not written and should remain zero at all times)
+Eps = np.zeros((nelem,nip,3,3))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # nodal displacements
+ ue = disp[nodes,:]
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ dNdxe = dNx[e,q,:,:]
+
+ # displacement gradient
+ gradu = np.einsum('mi,mj->ij', dNdxe, ue)
+
+ # compute strain tensor, and store per integration point
+ # (use plain strain to convert 2-d to 3-d tensor)
+ Eps[e,q,:2,:2] = .5 * ( gradu + gradu.T )
+
+# constitutive response: strain tensor -> stress tensor (per integration point)
+Sig = ddot42(C4, Eps)
+
+# ==================================================================================================
+# internal force from stress
+# ==================================================================================================
+
+# allocate internal force
+fint = np.zeros((ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
- Jinv = np.linalg.inv(J)
- Jdet = np.linalg.det(J)
+ # allocate element internal force
+ fe = np.zeros((nne*ndim))
- # -- shape function gradients (w.r.t. the global coordinates)
- dNdx = np.zeros((nne,nd))
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
- for m in range(nne):
- for i in range(nd):
- for j in range(nd):
- dNdx[m,i] += Jinv[i,j] * dNdxi[m,j]
+ # alias integration point values
+ # (plane strain: stress in z-direction irrelevant)
+ sig = Sig[e,q,:2,:2]
+ dNdxe = dNx[e,q,:,:]
+ dV = vol[e,q]
- # -- assemble to element stiffness matrix
- for m in range(nne):
- for n in range(nne):
- for i in range(nd):
- for j in range(nd):
- for k in range(nd):
- for l in range(nd):
- Ke[m*nd+j, n*nd+k] += dNdx[m,i] * C4[i,j,k,l] * dNdx[n,l] * w * Jdet * thick
+ # assemble to element internal force
+ # fe[m*nd+j] += dNdx[m,i] * sig[i,j] * dV
+ fe += ( np.einsum('mi,ij->mj', dNdxe, sig) * dV ).reshape(nne*ndim)
- # - assemble to global stiffness matrix
- iie = dofs[e,:].ravel()
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
+ fint[iie] += fe
+# ==================================================================================================
+# stiffness matrix from 'tangent'
+# ==================================================================================================
+
+# allocate stiffness matrix
+K = np.zeros((ndof, ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element stiffness matrix
+ Ke = np.zeros((nne*ndim, nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ c4 = C4 [e,q,:2,:2,:2,:2]
+ dNdxe = dNx[e,q,:,:]
+ dV = vol[e,q]
+
+ # assemble to element stiffness matrix
+ # Ke[m*nd+j, n*nd+k] += dNdx[m,i] * C4[i,j,k,l] * dNdx[n,l] * dV
+ Ke += ( np.einsum('mi,ijkl,nl->mjnk', dNdxe, c4, dNdxe) * dV ).reshape(nne*ndim, nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
K[np.ix_(iie,iie)] += Ke
# ==================================================================================================
# partition and solve
# ==================================================================================================
-# zero-initialize displacements and forces
-f = np.zeros((ndof))
-u = np.zeros((ndof))
+# prescribed external force: zero on all free DOFs
+# (other DOFs ignored in solution, the reaction forces on the DOFs are computed below)
+fext = np.zeros((ndof))
+
+# DOF-partitioning: ['u'nknown, 'p'rescribed]
+# - prescribed
+iip = np.hstack((
+ dofs[nodesBottom,1],
+ dofs[nodesLeft ,0],
+ dofs[nodesRight ,0],
+))
+# - unknown
+iiu = np.setdiff1d(dofs.ravel(), iip)
# fixed displacements
-# - zero-initialize
-iip = np.empty((0),dtype='int' )
-up = np.empty((0),dtype='float')
-# - 'y = 0' : 'u_y = 0'
-idof = dofs[inode[0,:], 1]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.0 * np.ones(idof.shape) ))
-# - 'x = 0' : 'u_x = 0'
-idof = dofs[inode[:,0], 0]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.0 * np.ones(idof.shape) ))
-# - 'x = 1' : 'u_x = 0.1'
-idof = dofs[inode[:,-1], 0]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.1 * np.ones(idof.shape) ))
-
-# free displacements
-iiu = np.setdiff1d(dofs.ravel(), iip)
+up = np.hstack((
+ 0.0 * np.ones((len(nodesBottom))),
+ 0.0 * np.ones((len(nodesLeft ))),
+ 0.1 * np.ones((len(nodesRight ))),
+))
+
+# residual force
+r = fext - fint
# partition
# - stiffness matrix
-Kuu = K[np.ix_(iiu, iiu)]
-Kup = K[np.ix_(iiu, iip)]
-Kpu = K[np.ix_(iip, iiu)]
-Kpp = K[np.ix_(iip, iip)]
-# - external force
-fu = f[iiu]
+Kuu = K[np.ix_(iiu, iiu)]
+Kup = K[np.ix_(iiu, iip)]
+Kpu = K[np.ix_(iip, iiu)]
+Kpp = K[np.ix_(iip, iip)]
+# - residual force
+ru = r[iiu]
-# solve
-uu = np.linalg.solve(Kuu, fu - Kup.dot(up))
-fp = Kpu.dot(uu) + Kpp.dot(up)
+# solve for unknown displacement DOFs
+uu = np.linalg.solve(Kuu, ru - Kup.dot(up))
# assemble
+u = np.empty((ndof))
u[iiu] = uu
u[iip] = up
-f[iip] = fp
-# convert to nodal displacement and forces
+# convert to nodal displacements
disp = u[dofs]
-fext = f[dofs]
+
+# ==================================================================================================
+# strain from nodal displacements, stress from constitutive response
+# ==================================================================================================
+
+# zero-initialise strain tensor per integration point
+# (plain strain -> all z-components are not written and should remain zero at all times)
+Eps = np.zeros((nelem,nip,3,3))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # nodal displacements
+ ue = disp[nodes,:]
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ dNdxe = dNx[e,q,:,:]
+
+ # displacement gradient
+ gradu = np.einsum('mi,mj->ij', dNdxe, ue)
+
+ # compute strain tensor, and store per integration point
+ # (use plain strain to convert 2-d to 3-d tensor)
+ Eps[e,q,:2,:2] = .5 * ( gradu + gradu.T )
+
+# constitutive response: strain tensor -> stress tensor (per integration point)
+Sig = ddot42(C4, Eps)
+
+# ==================================================================================================
+# internal force from stress, reaction (external) force from internal force on fixed nodes
+# ==================================================================================================
+
+# allocate internal force
+fint = np.zeros((ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element internal force
+ fe = np.zeros((nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ # (plane strain: stress in z-direction irrelevant)
+ sig = Sig[e,q,:2,:2]
+ dNdxe = dNx[e,q,:,:]
+ dV = vol[e,q]
+
+ # assemble to element internal force
+ # fe[m*nd+j] += dNdx[m,i] * sig[i,j] * dV
+ fe += ( np.einsum('mi,ij->mj', dNdxe, sig) * dV ).reshape(nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
+ fint[iie] += fe
+
+# reaction force
+fext[iip] = fint[iip]
# ==================================================================================================
# plot
# ==================================================================================================
import matplotlib.pyplot as plt
-fig, axes = plt.subplots(ncols=2, figsize=(16,8))
+fig, axes = plt.subplots(ncols=2, figsize=(16,8))
+
+# reconstruct external force as nodal vectors
+fext = fext[dofs]
+# plot original nodal positions + displacement field as arrows
axes[0].scatter(coor[:,0], coor[:,1])
axes[0].quiver (coor[:,0], coor[:,1], disp[:,0], disp[:,1])
+# plot new nodal positions + external (reaction) force field as arrows
axes[1].scatter((coor+disp)[:,0], (coor+disp)[:,1])
axes[1].quiver ((coor+disp)[:,0], (coor+disp)[:,1], fext[:,0], fext[:,1], scale=.1)
+# fix axes limits
for axis in axes:
axis.set_xlim([-0.4, 1.5])
axis.set_ylim([-0.4, 1.5])
plt.show()
diff --git a/docs/theory/examples/cartesian_linear-elasticity_B-matrix.py b/docs/theory/examples/cartesian_linear-elasticity_B-matrix.py
index 2fb964e..fb7e9dd 100644
--- a/docs/theory/examples/cartesian_linear-elasticity_B-matrix.py
+++ b/docs/theory/examples/cartesian_linear-elasticity_B-matrix.py
@@ -1,308 +1,392 @@
+# ==================================================================================================
+#
+# (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
+#
+# ==================================================================================================
import numpy as np
np.set_printoptions(linewidth=200)
# ==================================================================================================
# tensor products
# ==================================================================================================
-def ddot43(A,B):
-
- nd = A.shape[0]
-
- C = np.zeros((nd,nd,nd))
-
- for i in range(nd):
- for j in range(nd):
- for k in range(nd):
- for l in range(nd):
- for m in range(nd):
- C[i,j,m] += A[i,j,k,l] * B[l,k,m]
-
- return C
-
-# --------------------------------------------------------------------------------------------------
-
-def ddot33(A,B):
-
- nd = A.shape[0]
-
- C = np.zeros((nd,nd))
-
- for i in range(nd):
- for j in range(nd):
- for k in range(nd):
- for l in range(nd):
- C[i,l] += A[i,j,k] * B[k,j,l]
-
- return C
-
-# --------------------------------------------------------------------------------------------------
-
-def transpose3(A):
-
- nd = A.shape[0]
-
- C = np.zeros((nd,nd,nd))
-
- for i in range(nd):
- for j in range(nd):
- for k in range(nd):
- C[k,j,i] = A[i,j,k]
-
- return C
-
-# ==================================================================================================
-# identity tensors
-# ==================================================================================================
-
-II = np.zeros((3,3,3,3))
-I4 = np.zeros((3,3,3,3))
-I4rt = np.zeros((3,3,3,3))
-Is = np.zeros((3,3,3,3))
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == j and k == l ):
- II[i,j,k,l] = 1.
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == l and j == k ):
- I4[i,j,k,l] = 1.
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == k and j == l ):
- I4rt[i,j,k,l] = 1.
-
-I4s = ( I4 + I4rt ) / 2.
-I4d = ( I4s - II/3. )
-
-# ==================================================================================================
-# elasticity tensor
-# ==================================================================================================
-
-# bulk and shear modulus
-K = 0.8333333333333333
-G = 0.3846153846153846
-
-# elasticity tensor (3d)
-C4 = K * II + 2. * G * I4d
+ddot22 = lambda A2,B2: np.einsum('...ij ,...ji->... ',A2,B2)
+ddot42 = lambda A4,B2: np.einsum('...ijkl,...lk->...ij ',A4,B2)
+dyad22 = lambda A2,B2: np.einsum('...ij ,...kl->...ijkl',A2,B2)
# ==================================================================================================
# mesh definition
# ==================================================================================================
# number of elements
-nx = 10
-ny = 10
+nx = 20
+ny = 20
# mesh dimensions
nelem = nx * ny # number of elements
nnode = (nx+1) * (ny+1) # number of nodes
nne = 4 # number of nodes per element
-nd = 2 # number of dimensions
-ndof = nnode * nd # total number of degrees of freedom
+ndim = 2 # number of dimensions
+ndof = nnode * ndim # total number of degrees of freedom
# out-of-plane thickness
thick = 1.
-# coordinates and connectivity: zero-initialize
-coor = np.zeros((nnode,nd ), dtype='float')
-conn = np.zeros((nelem,nne), dtype='int' )
+# zero-initialise coordinates, displacements, and connectivity
+coor = np.zeros((nnode,ndim), dtype='float')
+disp = np.zeros((nnode,ndim), dtype='float')
+conn = np.zeros((nelem,nne ), dtype='int' )
-# coordinates: set
+# coordinates
# - grid of points
x,y = np.meshgrid(np.linspace(0,1,nx+1), np.linspace(0,1,ny+1))
# - store from grid of points
coor[:,0] = x.ravel()
coor[:,1] = y.ravel()
-# connectivity: set
+# connectivity
# - grid of node numbers
inode = np.arange(nnode).reshape(ny+1, nx+1)
# - store from grid of node numbers
conn[:,0] = inode[ :-1, :-1].ravel()
conn[:,1] = inode[ :-1,1: ].ravel()
conn[:,2] = inode[1: ,1: ].ravel()
conn[:,3] = inode[1: , :-1].ravel()
+# - node sets
+nodesLeft = inode[ :, 0]
+nodesRight = inode[ :,-1]
+nodesBottom = inode[ 0, :]
+nodesTop = inode[-1, :]
# DOF-numbers per node
-dofs = np.arange(nnode*nd).reshape(nnode,nd)
+dofs = np.arange(nnode*ndim).reshape(nnode,ndim)
# ==================================================================================================
# quadrature definition
# ==================================================================================================
# integration point coordinates (local element coordinates)
Xi = np.array([
[-1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), +1./np.sqrt(3.)],
[-1./np.sqrt(3.), +1./np.sqrt(3.)],
])
# integration point weights
W = np.array([
[1.],
[1.],
[1.],
[1.],
])
+# number of integration points
+nip = 4
+
# ==================================================================================================
-# assemble stiffness matrix
+# B-matrix at each integration point
# ==================================================================================================
-# allocate matrix
-K = np.zeros((ndof, ndof))
+# allocate
+B = np.empty((nelem,nip,nne,3,3,3))
+vol = np.empty((nelem,nip))
# loop over nodes
-for e in conn:
+for e, nodes in enumerate(conn):
- # - nodal coordinates
- xe = coor[e,:]
+ # nodal coordinates
+ xe = coor[nodes,:]
- # - allocate element stiffness matrix
- Ke = np.zeros((nne*nd, nne*nd))
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
- # - numerical quadrature
- for w, xi in zip(W, Xi):
-
- # -- shape function gradients (w.r.t. the local element coordinates)
+ # shape function gradients (w.r.t. the local element coordinates)
dNdxi = np.array([
[-.25*(1.-xi[1]), -.25*(1.-xi[0])],
[+.25*(1.-xi[1]), -.25*(1.+xi[0])],
[+.25*(1.+xi[1]), +.25*(1.+xi[0])],
[-.25*(1.+xi[1]), +.25*(1.-xi[0])],
])
- # -- Jacobian
- J = np.zeros((nd, nd))
+ # Jacobian
+ Je = np.einsum('mi,mj->ij', dNdxi, xe)
+ invJe = np.linalg.inv(Je)
+ detJe = np.linalg.det(Je)
+
+ # shape function gradients (w.r.t. the global coordinates)
+ dNdxe = np.einsum('ij,mj->mi', invJe, dNdxi)
+ # compute B-matrix and integration-point volume and store for later use
for m in range(nne):
- for i in range(nd):
- for j in range(nd):
- J[i,j] += dNdxi[m,i] * xe[m,j]
- Jinv = np.linalg.inv(J)
- Jdet = np.linalg.det(J)
+ Be = np.zeros((3,3,3))
- # -- shape function gradients (w.r.t. the global coordinates)
- dNdx = np.zeros((nne,nd))
+ Be[0,0,0] = dNdxe[m,0]
+ Be[0,1,1] = dNdxe[m,0]
+ Be[1,0,0] = dNdxe[m,1]
+ Be[1,1,1] = dNdxe[m,1]
- for m in range(nne):
- for i in range(nd):
- for j in range(nd):
- dNdx[m,i] += Jinv[i,j] * dNdxi[m,j]
+ B[e,q,m,:,:,:] = Be
- # -- assemble to element stiffness matrix
- for m in range(nne):
+ vol[e,q] = w * detJe * thick
+
+# ==================================================================================================
+# stiffness tensor at each integration point (provides constitutive response and 'tangent')
+# ==================================================================================================
+
+# identity tensors (per integration point)
+i = np.eye(3)
+I = np.einsum('ij ,...' , i ,np.ones([nelem,nip]))
+I4 = np.einsum('ijkl,...->...ijkl',np.einsum('il,jk',i,i),np.ones([nelem,nip]))
+I4rt = np.einsum('ijkl,...->...ijkl',np.einsum('ik,jl',i,i),np.ones([nelem,nip]))
+I4s = (I4+I4rt)/2.
+II = dyad22(I,I)
+I4d = I4s-II/3.
- Bm = np.zeros((3,3,3))
+# bulk and shear modulus (homogeneous)
+K = 0.8333333333333333 * np.ones((nelem, nip))
+G = 0.3846153846153846 * np.ones((nelem, nip))
- Bm[0,0,0] = dNdx[m,0]
- Bm[0,1,1] = dNdx[m,0]
+# elasticity tensor (per integration point)
+C4 = np.einsum('eq,eqijkl->eqijkl', K, II) + 2. * np.einsum('eq,eqijkl->eqijkl', G, I4d)
- Bm[1,0,0] = dNdx[m,1]
- Bm[1,1,1] = dNdx[m,1]
+# ==================================================================================================
+# strain from nodal displacements, stress from constitutive response
+# ==================================================================================================
- for n in range(nne):
+# allocate strain tensor per integration point
+Eps = np.empty((nelem,nip,3,3))
- Bn = np.zeros((3,3,3))
+# loop over nodes
+for e, nodes in enumerate(conn):
- Bn[0,0,0] = dNdx[n,0]
- Bn[0,1,1] = dNdx[n,0]
+ # nodal displacements in 3-d
+ # u_z = 0
+ ue = np.zeros((nne,3))
+ ue[:,0] = disp[nodes,0]
+ ue[:,1] = disp[nodes,1]
- Bn[1,0,0] = dNdx[n,1]
- Bn[1,1,1] = dNdx[n,1]
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
- Kmn = ddot33(transpose3(Bm),ddot43(C4,Bn))
+ # alias integration point values
+ Be = B[e,q,:,:,:,:]
- iim = np.array([ m*nd+0, m*nd+1 ])
- iin = np.array([ n*nd+0, n*nd+1 ])
+ # displacement gradient
+ gradu = np.einsum('mijk,mk->ij', Be, ue)
- Ke[np.ix_(iim,iin)] += Kmn[np.ix_([0,1], [0,1])] * w * Jdet * thick
+ # compute strain tensor, and store per integration point
+ Eps[e,q] = .5 * ( gradu + gradu.T )
- # - assemble to global stiffness matrix
- iie = dofs[e,:].ravel()
+# constitutive response: strain tensor -> stress tensor (per integration point)
+Sig = ddot42(C4, Eps)
+# ==================================================================================================
+# internal force from stress
+# ==================================================================================================
+
+# allocate internal force
+fint = np.zeros((ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element internal force
+ fe = np.zeros((nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ sig = Sig[e,q,:,:]
+ Be = B [e,q,:,:,:,:]
+ dV = vol[e,q]
+
+ # assemble to element internal force
+ # (plane strain: force in z-direction irrelevant)
+ # Bm = B[e,q,m,:,:,:]
+ # fm = ddot32(transpose3(Bm), sig) * dV
+ # fe[[m*ndim+0, m*ndim+1]] = [fm[m,2], fm[m,0]]
+ fq = np.einsum('mijk,ij->mk', Be, sig) * dV
+ fe += fq[:,:2].reshape(nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
+ fint[iie] += fe
+
+# ==================================================================================================
+# stiffness matrix from 'tangent'
+# ==================================================================================================
+
+# allocate stiffness matrix
+K = np.zeros((ndof, ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element stiffness matrix
+ Ke = np.zeros((nne*ndim, nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ c4 = C4 [e,q,:,:,:,:]
+ Be = B [e,q,:,:,:,:]
+ dV = vol[e,q]
+
+ # assemble to element stiffness matrix
+ # Bm = B[e,q,m,:,:,:]
+ # Bn = B[e,q,n,:,:,:]
+ # Kmn = ddot33(transpose3(Bm),ddot43(C4,Bn)) * dV
+ # Ke[[m*ndim+0,m*ndim+1], [n*ndim+0,n*ndim+1]] += Kmn[[2,0], [2,0]]
+ Kq = np.einsum('mabc,abde,nedf->mcnf', Be, c4, Be) * dV
+ Ke += Kq[:,:2,:,:2].reshape(nne*ndim, nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
K[np.ix_(iie,iie)] += Ke
# ==================================================================================================
# partition and solve
# ==================================================================================================
-# zero-initialize displacements and forces
-f = np.zeros((ndof))
-u = np.zeros((ndof))
+# prescribed external force: zero on all free DOFs
+# (other DOFs ignored in solution, the reaction forces on the DOFs are computed below)
+fext = np.zeros((ndof))
+
+# DOF-partitioning: ['u'nknown, 'p'rescribed]
+# - prescribed
+iip = np.hstack((
+ dofs[nodesBottom,1],
+ dofs[nodesLeft ,0],
+ dofs[nodesRight ,0],
+))
+# - unknown
+iiu = np.setdiff1d(dofs.ravel(), iip)
# fixed displacements
-# - zero-initialize
-iip = np.empty((0),dtype='int' )
-up = np.empty((0),dtype='float')
-# - 'y = 0' : 'u_y = 0'
-idof = dofs[inode[0,:], 1]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.0 * np.ones(idof.shape) ))
-# - 'x = 0' : 'u_x = 0'
-idof = dofs[inode[:,0], 0]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.0 * np.ones(idof.shape) ))
-# - 'x = 1' : 'u_x = 0.1'
-idof = dofs[inode[:,-1], 0]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.1 * np.ones(idof.shape) ))
-
-# free displacements
-iiu = np.setdiff1d(dofs.ravel(), iip)
+up = np.hstack((
+ 0.0 * np.ones((len(nodesBottom))),
+ 0.0 * np.ones((len(nodesLeft ))),
+ 0.1 * np.ones((len(nodesRight ))),
+))
+
+# residual force
+r = fext - fint
# partition
# - stiffness matrix
-Kuu = K[np.ix_(iiu, iiu)]
-Kup = K[np.ix_(iiu, iip)]
-Kpu = K[np.ix_(iip, iiu)]
-Kpp = K[np.ix_(iip, iip)]
-# - external force
-fu = f[iiu]
+Kuu = K[np.ix_(iiu, iiu)]
+Kup = K[np.ix_(iiu, iip)]
+Kpu = K[np.ix_(iip, iiu)]
+Kpp = K[np.ix_(iip, iip)]
+# - residual force
+ru = r[iiu]
-# solve
-uu = np.linalg.solve(Kuu, fu - Kup.dot(up))
-fp = Kpu.dot(uu) + Kpp.dot(up)
+# solve for unknown displacement DOFs
+uu = np.linalg.solve(Kuu, ru - Kup.dot(up))
# assemble
+u = np.empty((ndof))
u[iiu] = uu
u[iip] = up
-f[iip] = fp
-# convert to nodal displacement and forces
+# convert to nodal displacements
disp = u[dofs]
-fext = f[dofs]
+
+# ==================================================================================================
+# strain from nodal displacements, stress from constitutive response
+# ==================================================================================================
+
+# allocate strain tensor per integration point
+Eps = np.empty((nelem,nip,3,3))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # nodal displacements in 3-d
+ # u_z = 0
+ ue = np.zeros((nne,3))
+ ue[:,0] = disp[nodes,0]
+ ue[:,1] = disp[nodes,1]
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ Be = B[e,q,:,:,:,:]
+
+ # displacement gradient
+ gradu = np.einsum('mijk,mk->ij', Be, ue)
+
+ # compute strain tensor, and store per integration point
+ Eps[e,q] = .5 * ( gradu + gradu.T )
+
+# constitutive response: strain tensor -> stress tensor (per integration point)
+Sig = ddot42(C4, Eps)
+
+# ==================================================================================================
+# internal force from stress, reaction (external) force from internal force on fixed nodes
+# ==================================================================================================
+
+# allocate internal force
+fint = np.zeros((ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element internal force
+ fe = np.zeros((nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ sig = Sig[e,q,:,:]
+ Be = B [e,q,:,:,:,:]
+ dV = vol[e,q]
+
+ # assemble to element internal force
+ # (plane strain: force in z-direction irrelevant)
+ # Bm = B[e,q,m,:,:,:]
+ # fm = ddot32(transpose3(Bm), sig) * dV
+ # fe[[m*ndim+0, m*ndim+1]] = [fm[m,2], fm[m,0]]
+ fq = np.einsum('mijk,ij->mk', Be, sig) * dV
+ fe += fq[:,:2].reshape(nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
+ fint[iie] += fe
+
+# reaction force
+fext[iip] = fint[iip]
# ==================================================================================================
# plot
# ==================================================================================================
import matplotlib.pyplot as plt
-fig, axes = plt.subplots(ncols=2, figsize=(16,8))
+fig, axes = plt.subplots(ncols=2, figsize=(16,8))
+
+# reconstruct external force as nodal vectors
+fext = fext[dofs]
+# plot original nodal positions + displacement field as arrows
axes[0].scatter(coor[:,0], coor[:,1])
axes[0].quiver (coor[:,0], coor[:,1], disp[:,0], disp[:,1])
+# plot new nodal positions + external (reaction) force field as arrows
axes[1].scatter((coor+disp)[:,0], (coor+disp)[:,1])
axes[1].quiver ((coor+disp)[:,0], (coor+disp)[:,1], fext[:,0], fext[:,1], scale=.1)
+# fix axes limits
for axis in axes:
axis.set_xlim([-0.4, 1.5])
axis.set_ylim([-0.4, 1.5])
plt.show()
diff --git a/docs/theory/examples/cartesian_linear-elasticity_periodic-boundary-conditions.py b/docs/theory/examples/cartesian_linear-elasticity_periodic-boundary-conditions.py
new file mode 100644
index 0000000..82ed560
--- /dev/null
+++ b/docs/theory/examples/cartesian_linear-elasticity_periodic-boundary-conditions.py
@@ -0,0 +1,457 @@
+# ==================================================================================================
+#
+# (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
+#
+# ==================================================================================================
+
+import numpy as np
+
+np.set_printoptions(linewidth=200)
+
+# ==================================================================================================
+# tensor products
+# ==================================================================================================
+
+ddot22 = lambda A2,B2: np.einsum('...ij ,...ji->... ',A2,B2)
+ddot42 = lambda A4,B2: np.einsum('...ijkl,...lk->...ij ',A4,B2)
+dyad22 = lambda A2,B2: np.einsum('...ij ,...kl->...ijkl',A2,B2)
+
+# ==================================================================================================
+# mesh definition
+# ==================================================================================================
+
+# number of elements
+nx = 20
+ny = 20
+
+# mesh dimensions
+nelem = nx * ny # number of elements
+nnode = (nx+1) * (ny+1) # number of nodes
+nne = 4 # number of nodes per element
+ndim = 2 # number of dimensions
+ndof = nnode * ndim # total number of degrees of freedom
+
+# out-of-plane thickness
+thick = 1.
+
+# zero-initialise coordinates, displacements, and connectivity
+coor = np.zeros((nnode,ndim), dtype='float')
+disp = np.zeros((nnode,ndim), dtype='float')
+conn = np.zeros((nelem,nne ), dtype='int' )
+
+# coordinates
+# - grid of points
+x,y = np.meshgrid(np.linspace(0,1,nx+1), np.linspace(0,1,ny+1))
+# - store from grid of points
+coor[:,0] = x.ravel()
+coor[:,1] = y.ravel()
+
+# connectivity
+# - grid of node numbers
+inode = np.arange(nnode).reshape(ny+1, nx+1)
+# - store from grid of node numbers
+conn[:,0] = inode[ :-1, :-1].ravel()
+conn[:,1] = inode[ :-1,1: ].ravel()
+conn[:,2] = inode[1: ,1: ].ravel()
+conn[:,3] = inode[1: , :-1].ravel()
+# - node sets
+nodesLeft = inode[ :, 0]
+nodesRight = inode[ :,-1]
+nodesBottom = inode[ 0, :]
+nodesTop = inode[-1, :]
+
+# DOF-numbers per node
+dofs = np.arange(nnode*ndim).reshape(nnode,ndim)
+
+# ==================================================================================================
+# periodicity
+# ==================================================================================================
+
+# add virtual nodes
+# - node set
+nodesVirtual = nnode + np.arange(2)
+# - update size
+nnode += 2
+ndof += 4
+# - add coordinates (position is completely arbitrary)
+coor = np.vstack(( coor, np.array([[0.0, 1.1], [0.1, 1.1]]) ))
+# - add DOF numbers
+dofs = np.vstack(( dofs, np.max(dofs) + 1 + np.arange(4).reshape(2,2) ))
+
+# tyings (dependent, independent)
+tyings = np.vstack((
+ np.hstack(( nodesRight[1:-1].reshape(-1,1) , nodesLeft [1:-1].reshape(-1,1) )),
+ np.hstack(( nodesTop [1:-1].reshape(-1,1) , nodesBottom[1:-1].reshape(-1,1) )),
+ np.hstack(( nodesRight[ 0].reshape(-1,1) , nodesBottom[ 0].reshape(-1,1) )),
+ np.hstack(( nodesRight[ -1].reshape(-1,1) , nodesBottom[ 0].reshape(-1,1) )),
+ np.hstack(( nodesLeft [ -1].reshape(-1,1) , nodesBottom[ 0].reshape(-1,1) )),
+))
+
+# DOF-sets
+# - dependent DOFs
+iid = dofs[tyings[:,0],:].ravel()
+# - prescribed DOfs
+iip = np.hstack((
+ dofs[nodesBottom[0],:].ravel(),
+ dofs[nodesVirtual ,:].ravel(),
+))
+# - independent DOFs
+iii = np.setdiff1d(dofs.ravel(), iid)
+# - unknown DOFs
+iiu = np.setdiff1d(iii, iip)
+
+# renumber
+# - list
+renum = np.zeros(dofs.size, dtype='int')
+renum[iiu] = np.arange(iiu.size)
+renum[iip] = np.arange(iip.size) + iiu.size
+renum[iid] = np.arange(iid.size) + iiu.size + iip.size
+# - apply
+iiu = renum[iiu ]
+iip = renum[iip ]
+dofs = renum[dofs]
+# - define auxiliary DOF sets
+iii = np.hstack((iiu, iip))
+iid = np.arange(iid.size) + iii.size
+
+# nodal tyings
+# - zero-initialize
+Cdi = np.zeros((iid.size, iii.size))
+# - tie
+iie = 2 * np.arange(tyings.shape[0])
+Cdi[iie ,dofs[tyings[:,1],0]] = 1.
+Cdi[iie+1,dofs[tyings[:,1],1]] = 1.
+Cdi[iie ,dofs[nodesVirtual[0],0]] = coor[tyings[:,0],0] - coor[tyings[:,1],0] # (F-I)_xx
+Cdi[iie ,dofs[nodesVirtual[0],1]] = coor[tyings[:,0],1] - coor[tyings[:,1],1] # (F-I)_xy
+Cdi[iie+1,dofs[nodesVirtual[1],0]] = coor[tyings[:,0],0] - coor[tyings[:,1],0] # (F-I)_yx
+Cdi[iie+1,dofs[nodesVirtual[1],1]] = coor[tyings[:,0],1] - coor[tyings[:,1],1] # (F-I)_yy
+
+# ==================================================================================================
+# quadrature definition
+# ==================================================================================================
+
+# integration point coordinates (local element coordinates)
+Xi = np.array([
+ [-1./np.sqrt(3.), -1./np.sqrt(3.)],
+ [+1./np.sqrt(3.), -1./np.sqrt(3.)],
+ [+1./np.sqrt(3.), +1./np.sqrt(3.)],
+ [-1./np.sqrt(3.), +1./np.sqrt(3.)],
+])
+
+# integration point weights
+W = np.array([
+ [1.],
+ [1.],
+ [1.],
+ [1.],
+])
+
+# number of integration points
+nip = 4
+
+# ==================================================================================================
+# gradient of the shape functions at each integration point
+# ==================================================================================================
+
+# allocate
+dNx = np.empty((nelem,nip,nne,ndim))
+vol = np.empty((nelem,nip))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # nodal coordinates
+ xe = coor[nodes,:]
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # shape function gradients (w.r.t. the local element coordinates)
+ dNdxi = np.array([
+ [-.25*(1.-xi[1]), -.25*(1.-xi[0])],
+ [+.25*(1.-xi[1]), -.25*(1.+xi[0])],
+ [+.25*(1.+xi[1]), +.25*(1.+xi[0])],
+ [-.25*(1.+xi[1]), +.25*(1.-xi[0])],
+ ])
+
+ # Jacobian
+ Je = np.einsum('mi,mj->ij', dNdxi, xe)
+ invJe = np.linalg.inv(Je)
+ detJe = np.linalg.det(Je)
+
+ # shape function gradients (w.r.t. the global coordinates)
+ dNdxe = np.einsum('ij,mj->mi', invJe, dNdxi)
+
+ # store for later use
+ dNx[e,q,:,:] = dNdxe
+ vol[e,q] = w * detJe * thick
+
+# ==================================================================================================
+# stiffness tensor at each integration point (provides constitutive response and 'tangent')
+# ==================================================================================================
+
+# identity tensors (per integration point)
+i = np.eye(3)
+I = np.einsum('ij ,...' , i ,np.ones([nelem,nip]))
+I4 = np.einsum('ijkl,...->...ijkl',np.einsum('il,jk',i,i),np.ones([nelem,nip]))
+I4rt = np.einsum('ijkl,...->...ijkl',np.einsum('ik,jl',i,i),np.ones([nelem,nip]))
+I4s = (I4+I4rt)/2.
+II = dyad22(I,I)
+I4d = I4s-II/3.
+
+# bulk and shear modulus (homogeneous)
+K = 0.8333333333333333 * np.ones((nelem, nip))
+G = 0.3846153846153846 * np.ones((nelem, nip))
+
+# elasticity tensor (per integration point)
+C4 = np.einsum('eq,eqijkl->eqijkl', K, II) + 2. * np.einsum('eq,eqijkl->eqijkl', G, I4d)
+
+# ==================================================================================================
+# strain from nodal displacements, stress from constitutive response
+# ==================================================================================================
+
+# zero-initialise strain tensor per integration point
+# (plain strain -> all z-components are not written and should remain zero at all times)
+Eps = np.zeros((nelem,nip,3,3))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # nodal displacements
+ ue = disp[nodes,:]
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ dNdxe = dNx[e,q,:,:]
+
+ # displacement gradient
+ gradu = np.einsum('mi,mj->ij', dNdxe, ue)
+
+ # compute strain tensor, and store per integration point
+ # (use plain strain to convert 2-d to 3-d tensor)
+ Eps[e,q,:2,:2] = .5 * ( gradu + gradu.T )
+
+# constitutive response: strain tensor -> stress tensor (per integration point)
+Sig = ddot42(C4, Eps)
+
+# ==================================================================================================
+# internal force from stress
+# ==================================================================================================
+
+# allocate internal force
+fint = np.zeros((ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element internal force
+ fe = np.zeros((nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ # (plane strain: stress in z-direction irrelevant)
+ sig = Sig[e,q,:2,:2]
+ dNdxe = dNx[e,q,:,:]
+ dV = vol[e,q]
+
+ # assemble to element internal force
+ # fe[m*nd+j] += dNdx[m,i] * sig[i,j] * dV
+ fe += ( np.einsum('mi,ij->mj', dNdxe, sig) * dV ).reshape(nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
+ fint[iie] += fe
+
+# ==================================================================================================
+# stiffness matrix from 'tangent'
+# ==================================================================================================
+
+# allocate stiffness matrix
+K = np.zeros((ndof, ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element stiffness matrix
+ Ke = np.zeros((nne*ndim, nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ c4 = C4 [e,q,:2,:2,:2,:2]
+ dNdxe = dNx[e,q,:,:]
+ dV = vol[e,q]
+
+ # assemble to element stiffness matrix
+ # Ke[m*nd+j, n*nd+k] += dNdx[m,i] * C4[i,j,k,l] * dNdx[n,l] * dV
+ Ke += ( np.einsum('mi,ijkl,nl->mjnk', dNdxe, c4, dNdxe) * dV ).reshape(nne*ndim, nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
+ K[np.ix_(iie,iie)] += Ke
+
+# ==================================================================================================
+# partition and solve
+# ==================================================================================================
+
+# prescribed external force: zero on all free DOFs
+# (other DOFs ignored in solution, the reaction forces on the DOFs are computed below)
+fext = np.zeros((ndof))
+
+# fixed displacements
+up = np.array([
+ 0.0, # suppress rigid-body motion in x
+ 0.0, # suppress rigid-body motion in x
+ 0.0, # (F-I)_xx
+ 0.1, # (F-I)_xy
+ 0.0, # (F-I)_yx
+ 0.0, # (F-I)_yy
+])
+
+# residual force
+r = fext - fint
+
+# partition in independent and dependent part
+# - stiffness matrix
+Kii = K[np.ix_(iii, iii)]
+Kid = K[np.ix_(iii, iid)]
+Kdi = K[np.ix_(iid, iii)]
+Kdd = K[np.ix_(iid, iid)]
+# - residual force
+ri = r[iii]
+rd = r[iid]
+
+# condense system
+Kii = Kii + np.dot(Kid, Cdi) + np.dot(Cdi.T, Kdi) + np.dot(Cdi.T, np.dot(Kdd, Cdi))
+ri = ri + np.dot(Cdi.T, rd)
+
+# partition
+# - stiffness matrix
+Kuu = Kii[np.ix_(iiu, iiu)]
+Kup = Kii[np.ix_(iiu, iip)]
+Kpu = Kii[np.ix_(iip, iiu)]
+Kpp = Kii[np.ix_(iip, iip)]
+# - residual force
+ru = ri[iiu]
+
+# solve for unknown displacement DOFs
+uu = np.linalg.solve(Kuu, ru - Kup.dot(up))
+
+# assemble
+ui = np.empty((iii.size))
+ui[iiu] = uu
+ui[iip] = up
+
+# reconstruct
+# - dependent displacement DOFs
+ud = np.dot(Cdi, ui)
+# - assemble
+u = np.empty((ndof))
+u[iii] = ui
+u[iid] = ud
+
+# convert to nodal displacements
+disp = u[dofs]
+
+# ==================================================================================================
+# strain from nodal displacements, stress from constitutive response
+# ==================================================================================================
+
+# zero-initialise strain tensor per integration point
+# (plain strain -> all z-components are not written and should remain zero at all times)
+Eps = np.zeros((nelem,nip,3,3))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # nodal displacements
+ ue = disp[nodes,:]
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ dNdxe = dNx[e,q,:,:]
+
+ # displacement gradient
+ gradu = np.einsum('mi,mj->ij', dNdxe, ue)
+
+ # compute strain tensor, and store per integration point
+ # (use plain strain to convert 2-d to 3-d tensor)
+ Eps[e,q,:2,:2] = .5 * ( gradu + gradu.T )
+
+# constitutive response: strain tensor -> stress tensor (per integration point)
+Sig = ddot42(C4, Eps)
+
+# ==================================================================================================
+# internal force from stress, reaction (external) force from internal force on fixed nodes
+# ==================================================================================================
+
+# allocate internal force
+fint = np.zeros((ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element internal force
+ fe = np.zeros((nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ # (plane strain: stress in z-direction irrelevant)
+ sig = Sig[e,q,:2,:2]
+ dNdxe = dNx[e,q,:,:]
+ dV = vol[e,q]
+
+ # assemble to element internal force
+ # fe[m*nd+j] += dNdx[m,i] * sig[i,j] * dV
+ fe += ( np.einsum('mi,ij->mj', dNdxe, sig) * dV ).reshape(nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
+ fint[iie] += fe
+
+# residual force
+# - all DOFs
+r = fext - fint
+# - partition [independent, dependent]
+ri = r[iii]
+rd = r[iid]
+# - condense system: apply tying relations
+ri = ri + np.dot(Cdi.T, rd)
+
+# reaction force
+fext[iip] = -ri[iip]
+
+# ==================================================================================================
+# plot
+# ==================================================================================================
+
+import matplotlib.pyplot as plt
+
+fig, axes = plt.subplots(ncols=2, figsize=(16,8))
+
+# reconstruct external force as nodal vectors (for all nodes)
+fext = fext[dofs]
+
+# plot original nodal positions + displacement field as arrows
+axes[0].scatter(coor[:,0], coor[:,1])
+axes[0].quiver (coor[:,0], coor[:,1], disp[:,0], disp[:,1])
+
+# plot new nodal positions + external (reaction) force field as arrows
+axes[1].scatter((coor+disp)[:,0], (coor+disp)[:,1])
+axes[1].quiver ((coor+disp)[:,0], (coor+disp)[:,1], fext[:,0], fext[:,1], scale=.1)
+
+# fix axes limits
+for axis in axes:
+ axis.set_xlim([-0.4, 1.5])
+ axis.set_ylim([-0.4, 1.5])
+
+plt.show()
diff --git a/docs/theory/examples/cartesian_linear-elasticity_sparse.py b/docs/theory/examples/cartesian_linear-elasticity_sparse.py
index 221df5e..670455e 100644
--- a/docs/theory/examples/cartesian_linear-elasticity_sparse.py
+++ b/docs/theory/examples/cartesian_linear-elasticity_sparse.py
@@ -1,241 +1,372 @@
+# ==================================================================================================
+#
+# (c - GPLv3) T.W.J. de Geus (Tom) | tom@geus.me | www.geus.me | github.com/tdegeus/GooseFEM
+#
+# ==================================================================================================
import numpy as np
from scipy.sparse import dok_matrix
from scipy.sparse.linalg import spsolve
np.set_printoptions(linewidth=200)
# ==================================================================================================
-# identity tensors
-# ==================================================================================================
-
-II = np.zeros((3,3,3,3))
-I4 = np.zeros((3,3,3,3))
-I4rt = np.zeros((3,3,3,3))
-Is = np.zeros((3,3,3,3))
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == j and k == l ):
- II[i,j,k,l] = 1.
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == l and j == k ):
- I4[i,j,k,l] = 1.
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == k and j == l ):
- I4rt[i,j,k,l] = 1.
-
-I4s = ( I4 + I4rt ) / 2.
-I4d = ( I4s - II/3. )
-
-# ==================================================================================================
-# elasticity tensor
+# tensor products
# ==================================================================================================
-# bulk and shear modulus
-K = 0.8333333333333333
-G = 0.3846153846153846
-
-# elasticity tensor (3d)
-C4 = K * II + 2. * G * I4d
+ddot22 = lambda A2,B2: np.einsum('...ij ,...ji->... ',A2,B2)
+ddot42 = lambda A4,B2: np.einsum('...ijkl,...lk->...ij ',A4,B2)
+dyad22 = lambda A2,B2: np.einsum('...ij ,...kl->...ijkl',A2,B2)
# ==================================================================================================
# mesh definition
# ==================================================================================================
# number of elements
nx = 20
ny = 20
# mesh dimensions
nelem = nx * ny # number of elements
nnode = (nx+1) * (ny+1) # number of nodes
nne = 4 # number of nodes per element
-nd = 2 # number of dimensions
-ndof = nnode * nd # total number of degrees of freedom
+ndim = 2 # number of dimensions
+ndof = nnode * ndim # total number of degrees of freedom
# out-of-plane thickness
thick = 1.
-# coordinates and connectivity: zero-initialize
-coor = np.zeros((nnode,nd ), dtype='float')
-conn = np.zeros((nelem,nne), dtype='int' )
+# zero-initialise coordinates, displacements, and connectivity
+coor = np.zeros((nnode,ndim), dtype='float')
+disp = np.zeros((nnode,ndim), dtype='float')
+conn = np.zeros((nelem,nne ), dtype='int' )
-# coordinates: set
+# coordinates
# - grid of points
x,y = np.meshgrid(np.linspace(0,1,nx+1), np.linspace(0,1,ny+1))
# - store from grid of points
coor[:,0] = x.ravel()
coor[:,1] = y.ravel()
-# connectivity: set
+# connectivity
# - grid of node numbers
inode = np.arange(nnode).reshape(ny+1, nx+1)
# - store from grid of node numbers
conn[:,0] = inode[ :-1, :-1].ravel()
conn[:,1] = inode[ :-1,1: ].ravel()
conn[:,2] = inode[1: ,1: ].ravel()
conn[:,3] = inode[1: , :-1].ravel()
+# - node sets
+nodesLeft = inode[ :, 0]
+nodesRight = inode[ :,-1]
+nodesBottom = inode[ 0, :]
+nodesTop = inode[-1, :]
# DOF-numbers per node
-dofs = np.arange(nnode*nd).reshape(nnode,nd)
+dofs = np.arange(nnode*ndim).reshape(nnode,ndim)
# ==================================================================================================
# quadrature definition
# ==================================================================================================
# integration point coordinates (local element coordinates)
Xi = np.array([
[-1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), -1./np.sqrt(3.)],
[+1./np.sqrt(3.), +1./np.sqrt(3.)],
[-1./np.sqrt(3.), +1./np.sqrt(3.)],
])
# integration point weights
W = np.array([
[1.],
[1.],
[1.],
[1.],
])
+# number of integration points
+nip = 4
+
# ==================================================================================================
-# assemble stiffness matrix
+# gradient of the shape functions at each integration point
# ==================================================================================================
-# allocate matrix
-K = dok_matrix((ndof, ndof), dtype=np.float)
+# allocate
+dNx = np.empty((nelem,nip,nne,ndim))
+vol = np.empty((nelem,nip))
# loop over nodes
-for e in conn:
+for e, nodes in enumerate(conn):
- # - nodal coordinates
- xe = coor[e,:]
+ # nodal coordinates
+ xe = coor[nodes,:]
- # - allocate element stiffness matrix
- Ke = np.zeros((nne*nd, nne*nd))
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
- # - numerical quadrature
- for w, xi in zip(W, Xi):
-
- # -- shape function gradients (w.r.t. the local element coordinates)
+ # shape function gradients (w.r.t. the local element coordinates)
dNdxi = np.array([
[-.25*(1.-xi[1]), -.25*(1.-xi[0])],
[+.25*(1.-xi[1]), -.25*(1.+xi[0])],
[+.25*(1.+xi[1]), +.25*(1.+xi[0])],
[-.25*(1.+xi[1]), +.25*(1.-xi[0])],
])
- # -- Jacobian
- J = np.zeros((nd, nd))
+ # Jacobian
+ Je = np.einsum('mi,mj->ij', dNdxi, xe)
+ invJe = np.linalg.inv(Je)
+ detJe = np.linalg.det(Je)
+
+ # shape function gradients (w.r.t. the global coordinates)
+ dNdxe = np.einsum('ij,mj->mi', invJe, dNdxi)
+
+ # store for later use
+ dNx[e,q,:,:] = dNdxe
+ vol[e,q] = w * detJe * thick
+
+# ==================================================================================================
+# stiffness tensor at each integration point (provides constitutive response and 'tangent')
+# ==================================================================================================
+
+# identity tensors (per integration point)
+i = np.eye(3)
+I = np.einsum('ij ,...' , i ,np.ones([nelem,nip]))
+I4 = np.einsum('ijkl,...->...ijkl',np.einsum('il,jk',i,i),np.ones([nelem,nip]))
+I4rt = np.einsum('ijkl,...->...ijkl',np.einsum('ik,jl',i,i),np.ones([nelem,nip]))
+I4s = (I4+I4rt)/2.
+II = dyad22(I,I)
+I4d = I4s-II/3.
- for m in range(nne):
- for i in range(nd):
- for j in range(nd):
- J[i,j] += dNdxi[m,i] * xe[m,j]
+# bulk and shear modulus (homogeneous)
+K = 0.8333333333333333 * np.ones((nelem, nip))
+G = 0.3846153846153846 * np.ones((nelem, nip))
+
+# elasticity tensor (per integration point)
+C4 = np.einsum('eq,eqijkl->eqijkl', K, II) + 2. * np.einsum('eq,eqijkl->eqijkl', G, I4d)
+
+# ==================================================================================================
+# strain from nodal displacements, stress from constitutive response
+# ==================================================================================================
+
+# zero-initialise strain tensor per integration point
+# (plain strain -> all z-components are not written and should remain zero at all times)
+Eps = np.zeros((nelem,nip,3,3))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # nodal displacements
+ ue = disp[nodes,:]
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ dNdxe = dNx[e,q,:,:]
+
+ # displacement gradient
+ gradu = np.einsum('mi,mj->ij', dNdxe, ue)
+
+ # compute strain tensor, and store per integration point
+ # (use plain strain to convert 2-d to 3-d tensor)
+ Eps[e,q,:2,:2] = .5 * ( gradu + gradu.T )
+
+# constitutive response: strain tensor -> stress tensor (per integration point)
+Sig = ddot42(C4, Eps)
+
+# ==================================================================================================
+# internal force from stress
+# ==================================================================================================
+
+# allocate internal force
+fint = np.zeros((ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
- Jinv = np.linalg.inv(J)
- Jdet = np.linalg.det(J)
+ # allocate element internal force
+ fe = np.zeros((nne*ndim))
- # -- shape function gradients (w.r.t. the global coordinates)
- dNdx = np.zeros((nne,nd))
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
- for m in range(nne):
- for i in range(nd):
- for j in range(nd):
- dNdx[m,i] += Jinv[i,j] * dNdxi[m,j]
+ # alias integration point values
+ # (plane strain: stress in z-direction irrelevant)
+ sig = Sig[e,q,:2,:2]
+ dNdxe = dNx[e,q,:,:]
+ dV = vol[e,q]
- # -- assemble to element stiffness matrix
- for m in range(nne):
- for n in range(nne):
- for i in range(nd):
- for j in range(nd):
- for k in range(nd):
- for l in range(nd):
- Ke[m*nd+j, n*nd+k] += dNdx[m,i] * C4[i,j,k,l] * dNdx[n,l] * w * Jdet * thick
+ # assemble to element internal force
+ # fe[m*nd+j] += dNdx[m,i] * sig[i,j] * dV
+ fe += ( np.einsum('mi,ij->mj', dNdxe, sig) * dV ).reshape(nne*ndim)
- # - assemble to global stiffness matrix
- iie = dofs[e,:].ravel()
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
+ fint[iie] += fe
+# ==================================================================================================
+# stiffness matrix from 'tangent'
+# ==================================================================================================
+
+# allocate stiffness matrix
+K = dok_matrix((ndof, ndof), dtype=np.float)
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element stiffness matrix
+ Ke = np.zeros((nne*ndim, nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ c4 = C4 [e,q,:2,:2,:2,:2]
+ dNdxe = dNx[e,q,:,:]
+ dV = vol[e,q]
+
+ # assemble to element stiffness matrix
+ # Ke[m*nd+j, n*nd+k] += dNdx[m,i] * C4[i,j,k,l] * dNdx[n,l] * dV
+ Ke += ( np.einsum('mi,ijkl,nl->mjnk', dNdxe, c4, dNdxe) * dV ).reshape(nne*ndim, nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
K[np.ix_(iie,iie)] += Ke
# ==================================================================================================
# partition and solve
# ==================================================================================================
-# zero-initialize displacements and forces
-f = np.zeros((ndof))
-u = np.zeros((ndof))
+# prescribed external force: zero on all free DOFs
+# (other DOFs ignored in solution, the reaction forces on the DOFs are computed below)
+fext = np.zeros((ndof))
+
+# DOF-partitioning: ['u'nknown, 'p'rescribed]
+# - prescribed
+iip = np.hstack((
+ dofs[nodesBottom,1],
+ dofs[nodesLeft ,0],
+ dofs[nodesRight ,0],
+))
+# - unknown
+iiu = np.setdiff1d(dofs.ravel(), iip)
# fixed displacements
-# - zero-initialize
-iip = np.empty((0),dtype='int' )
-up = np.empty((0),dtype='float')
-# - 'y = 0' : 'u_y = 0'
-idof = dofs[inode[0,:], 1]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.0 * np.ones(idof.shape) ))
-# - 'x = 0' : 'u_x = 0'
-idof = dofs[inode[:,0], 0]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.0 * np.ones(idof.shape) ))
-# - 'x = 1' : 'u_x = 0.1'
-idof = dofs[inode[:,-1], 0]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.1 * np.ones(idof.shape) ))
-
-# free displacements
-iiu = np.setdiff1d(dofs.ravel(), iip)
+up = np.hstack((
+ 0.0 * np.ones((len(nodesBottom))),
+ 0.0 * np.ones((len(nodesLeft ))),
+ 0.1 * np.ones((len(nodesRight ))),
+))
+
+# residual force
+r = fext - fint
# partition
# - stiffness matrix
Kuu = K.tocsr()[iiu,:].tocsc()[:,iiu]
Kup = K.tocsr()[iiu,:].tocsc()[:,iip]
Kpu = K.tocsr()[iip,:].tocsc()[:,iiu]
Kpp = K.tocsr()[iip,:].tocsc()[:,iip]
-# - external force
-fu = f[iiu]
+# - residual force
+ru = r[iiu]
-# solve
-uu = spsolve(Kuu, fu - Kup.dot(up))
-fp = Kpu.dot(uu) + Kpp.dot(up)
+# solve for unknown displacement DOFs
+uu = spsolve(Kuu, ru - Kup.dot(up))
# assemble
+u = np.empty((ndof))
u[iiu] = uu
u[iip] = up
-f[iip] = fp
-# convert to nodal displacement and forces
+# convert to nodal displacements
disp = u[dofs]
-fext = f[dofs]
+
+# ==================================================================================================
+# strain from nodal displacements, stress from constitutive response
+# ==================================================================================================
+
+# zero-initialise strain tensor per integration point
+# (plain strain -> all z-components are not written and should remain zero at all times)
+Eps = np.zeros((nelem,nip,3,3))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # nodal displacements
+ ue = disp[nodes,:]
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ dNdxe = dNx[e,q,:,:]
+
+ # displacement gradient
+ gradu = np.einsum('mi,mj->ij', dNdxe, ue)
+
+ # compute strain tensor, and store per integration point
+ # (use plain strain to convert 2-d to 3-d tensor)
+ Eps[e,q,:2,:2] = .5 * ( gradu + gradu.T )
+
+# constitutive response: strain tensor -> stress tensor (per integration point)
+Sig = ddot42(C4, Eps)
+
+# ==================================================================================================
+# internal force from stress, reaction (external) force from internal force on fixed nodes
+# ==================================================================================================
+
+# allocate internal force
+fint = np.zeros((ndof))
+
+# loop over nodes
+for e, nodes in enumerate(conn):
+
+ # allocate element internal force
+ fe = np.zeros((nne*ndim))
+
+ # loop over integration points
+ for q, (w, xi) in enumerate(zip(W, Xi)):
+
+ # alias integration point values
+ # (plane strain: stress in z-direction irrelevant)
+ sig = Sig[e,q,:2,:2]
+ dNdxe = dNx[e,q,:,:]
+ dV = vol[e,q]
+
+ # assemble to element internal force
+ # fe[m*nd+j] += dNdx[m,i] * sig[i,j] * dV
+ fe += ( np.einsum('mi,ij->mj', dNdxe, sig) * dV ).reshape(nne*ndim)
+
+ # assemble to global stiffness matrix
+ iie = dofs[nodes,:].ravel()
+ fint[iie] += fe
+
+# reaction force
+fext[iip] = fint[iip]
# ==================================================================================================
# plot
# ==================================================================================================
import matplotlib.pyplot as plt
fig, axes = plt.subplots(ncols=2, figsize=(16,8))
+# reconstruct external force as nodal vectors
+fext = fext[dofs]
+
+# plot original nodal positions + displacement field as arrows
axes[0].scatter(coor[:,0], coor[:,1])
axes[0].quiver (coor[:,0], coor[:,1], disp[:,0], disp[:,1])
+# plot new nodal positions + external (reaction) force field as arrows
axes[1].scatter((coor+disp)[:,0], (coor+disp)[:,1])
axes[1].quiver ((coor+disp)[:,0], (coor+disp)[:,1], fext[:,0], fext[:,1], scale=.1)
+# fix axes limits
for axis in axes:
axis.set_xlim([-0.4, 1.5])
axis.set_ylim([-0.4, 1.5])
plt.show()
diff --git a/docs/theory/examples/cartesian_linear-elasticity_sparse_strain.py b/docs/theory/examples/cartesian_linear-elasticity_sparse_strain.py
deleted file mode 100644
index 5827bc9..0000000
--- a/docs/theory/examples/cartesian_linear-elasticity_sparse_strain.py
+++ /dev/null
@@ -1,302 +0,0 @@
-
-import numpy as np
-from scipy.sparse import dok_matrix
-from scipy.sparse.linalg import spsolve
-
-np.set_printoptions(linewidth=200)
-
-# ==================================================================================================
-# identity tensors
-# ==================================================================================================
-
-II = np.zeros((3,3,3,3))
-I4 = np.zeros((3,3,3,3))
-I4rt = np.zeros((3,3,3,3))
-Is = np.zeros((3,3,3,3))
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == j and k == l ):
- II[i,j,k,l] = 1.
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == l and j == k ):
- I4[i,j,k,l] = 1.
-
-for i in range(3):
- for j in range(3):
- for k in range(3):
- for l in range(3):
- if ( i == k and j == l ):
- I4rt[i,j,k,l] = 1.
-
-I4s = ( I4 + I4rt ) / 2.
-I4d = ( I4s - II/3. )
-
-# ==================================================================================================
-# elasticity tensor
-# ==================================================================================================
-
-# bulk and shear modulus
-K = 0.8333333333333333
-G = 0.3846153846153846
-
-# elasticity tensor (3d)
-C4 = K * II + 2. * G * I4d
-
-# ==================================================================================================
-# mesh definition
-# ==================================================================================================
-
-# number of elements
-nx = 20
-ny = 20
-
-# mesh dimensions
-nelem = nx * ny # number of elements
-nnode = (nx+1) * (ny+1) # number of nodes
-nne = 4 # number of nodes per element
-nd = 2 # number of dimensions
-nip = 4 # number of integration points
-ndof = nnode * nd # total number of degrees of freedom
-
-# out-of-plane thickness
-thick = 1.
-
-# coordinates and connectivity: zero-initialize
-coor = np.zeros((nnode,nd ), dtype='float')
-conn = np.zeros((nelem,nne), dtype='int' )
-
-# coordinates: set
-# - grid of points
-x,y = np.meshgrid(np.linspace(0,1,nx+1), np.linspace(0,1,ny+1))
-# - store from grid of points
-coor[:,0] = x.ravel()
-coor[:,1] = y.ravel()
-
-# connectivity: set
-# - grid of node numbers
-inode = np.arange(nnode).reshape(ny+1, nx+1)
-# - store from grid of node numbers
-conn[:,0] = inode[ :-1, :-1].ravel()
-conn[:,1] = inode[ :-1,1: ].ravel()
-conn[:,2] = inode[1: ,1: ].ravel()
-conn[:,3] = inode[1: , :-1].ravel()
-
-# DOF-numbers per node
-dofs = np.arange(nnode*nd).reshape(nnode,nd)
-
-# ==================================================================================================
-# quadrature definition
-# ==================================================================================================
-
-# integration point coordinates (local element coordinates)
-Xi = np.array([
- [-1./np.sqrt(3.), -1./np.sqrt(3.)],
- [+1./np.sqrt(3.), -1./np.sqrt(3.)],
- [+1./np.sqrt(3.), +1./np.sqrt(3.)],
- [-1./np.sqrt(3.), +1./np.sqrt(3.)],
-])
-
-# integration point weights
-W = np.array([
- [1.],
- [1.],
- [1.],
- [1.],
-])
-
-# ==================================================================================================
-# assemble stiffness matrix
-# ==================================================================================================
-
-# allocate matrix
-K = dok_matrix((ndof, ndof), dtype=np.float)
-
-# loop over nodes
-for e in conn:
-
- # - nodal coordinates
- xe = coor[e,:]
-
- # - allocate element stiffness matrix
- Ke = np.zeros((nne*nd, nne*nd))
-
- # - numerical quadrature
- for w, xi in zip(W, Xi):
-
- # -- shape function gradients (w.r.t. the local element coordinates)
- dNdxi = np.array([
- [-.25*(1.-xi[1]), -.25*(1.-xi[0])],
- [+.25*(1.-xi[1]), -.25*(1.+xi[0])],
- [+.25*(1.+xi[1]), +.25*(1.+xi[0])],
- [-.25*(1.+xi[1]), +.25*(1.-xi[0])],
- ])
-
- # -- Jacobian
- J = np.zeros((nd, nd))
-
- for m in range(nne):
- for i in range(nd):
- for j in range(nd):
- J[i,j] += dNdxi[m,i] * xe[m,j]
-
- Jinv = np.linalg.inv(J)
- Jdet = np.linalg.det(J)
-
- # -- shape function gradients (w.r.t. the global coordinates)
- dNdx = np.zeros((nne,nd))
-
- for m in range(nne):
- for i in range(nd):
- for j in range(nd):
- dNdx[m,i] += Jinv[i,j] * dNdxi[m,j]
-
- # -- assemble to element stiffness matrix
- for m in range(nne):
- for n in range(nne):
- for i in range(nd):
- for j in range(nd):
- for k in range(nd):
- for l in range(nd):
- Ke[m*nd+j, n*nd+k] += dNdx[m,i] * C4[i,j,k,l] * dNdx[n,l] * w * Jdet * thick
-
- # - assemble to global stiffness matrix
- iie = dofs[e,:].ravel()
-
- K[np.ix_(iie,iie)] += Ke
-
-# ==================================================================================================
-# partition and solve
-# ==================================================================================================
-
-# zero-initialize displacements and forces
-f = np.zeros((ndof))
-u = np.zeros((ndof))
-
-# fixed displacements
-# - zero-initialize
-iip = np.empty((0),dtype='int' )
-up = np.empty((0),dtype='float')
-# - 'y = 0' : 'u_y = 0'
-idof = dofs[inode[0,:], 1]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.0 * np.ones(idof.shape) ))
-# - 'x = 0' : 'u_x = 0'
-idof = dofs[inode[:,0], 0]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.0 * np.ones(idof.shape) ))
-# - 'x = 1' : 'u_x = 0.1'
-idof = dofs[inode[:,-1], 0]
-iip = np.hstack(( iip, idof ))
-up = np.hstack(( up , 0.1 * np.ones(idof.shape) ))
-
-# free displacements
-iiu = np.setdiff1d(dofs.ravel(), iip)
-
-# partition
-# - stiffness matrix
-Kuu = K.tocsr()[iiu,:].tocsc()[:,iiu]
-Kup = K.tocsr()[iiu,:].tocsc()[:,iip]
-Kpu = K.tocsr()[iip,:].tocsc()[:,iiu]
-Kpp = K.tocsr()[iip,:].tocsc()[:,iip]
-# - external force
-fu = f[iiu]
-
-# solve
-uu = spsolve(Kuu, fu - Kup.dot(up))
-fp = Kpu.dot(uu) + Kpp.dot(up)
-
-# assemble
-u[iiu] = uu
-u[iip] = up
-f[iip] = fp
-
-# convert to nodal displacement and forces
-disp = u[dofs]
-fext = f[dofs]
-
-# ==================================================================================================
-# compute strain
-# ==================================================================================================
-
-# strain per integration point: allocate
-Eps = np.zeros((nelem,nip,nd,nd),dtype='float')
-
-# loop over nodes
-for ielem, e in enumerate(conn):
-
- # - nodal coordinates and displacements
- xe = coor[e,:]
- ue = disp[e,:]
-
- # - allocate element stiffness matrix
- Ke = np.zeros((nne*nd, nne*nd))
-
- # - numerical quadrature
- for k, (w, xi) in enumerate(zip(W, Xi)):
-
- # -- shape function gradients (w.r.t. the local element coordinates)
- dNdxi = np.array([
- [-.25*(1.-xi[1]), -.25*(1.-xi[0])],
- [+.25*(1.-xi[1]), -.25*(1.+xi[0])],
- [+.25*(1.+xi[1]), +.25*(1.+xi[0])],
- [-.25*(1.+xi[1]), +.25*(1.-xi[0])],
- ])
-
- # -- Jacobian
- J = np.zeros((nd, nd))
-
- for m in range(nne):
- for i in range(nd):
- for j in range(nd):
- J[i,j] += dNdxi[m,i] * xe[m,j]
-
- Jinv = np.linalg.inv(J)
- Jdet = np.linalg.det(J)
-
- # -- shape function gradients (w.r.t. the global coordinates)
- dNdx = np.zeros((nne,nd))
-
- for m in range(nne):
- for i in range(nd):
- for j in range(nd):
- dNdx[m,i] += Jinv[i,j] * dNdxi[m,j]
-
- # -- displacement gradient
- gradu = np.zeros((nd,nd))
-
- for m in range(nne):
- for i in range(nd):
- for j in range(nd):
- gradu[i,j] += dNdx[m,i] * ue[m,j]
-
- # -- strain
- eps = .5 * ( gradu + gradu.T )
-
- # -- assemble
- Eps[ielem,k,:,:] = eps
-
-# ==================================================================================================
-# plot
-# ==================================================================================================
-
-import matplotlib.pyplot as plt
-
-fig, axes = plt.subplots(ncols=2, figsize=(16,8))
-
-axes[0].scatter(coor[:,0], coor[:,1])
-axes[0].quiver (coor[:,0], coor[:,1], disp[:,0], disp[:,1])
-
-axes[1].imshow(Eps[:,:,0,0].reshape(ny*int(nip/nd), nx*int(nip/nd)), clim=(0,.2))
-
-for axis in axes:
- axis.set_xlim([-0.4, 1.5])
- axis.set_ylim([-0.4, 1.5])
-
-plt.show()
diff --git a/docs/theory/goose-article.cls b/docs/theory/goose-article.cls
index c3484e4..efe05b0 100644
--- a/docs/theory/goose-article.cls
+++ b/docs/theory/goose-article.cls
@@ -1,290 +1,321 @@
% ------------------------------------------------------------------------------
% @LaTeX-class-file{
% filename = "goose-article.cls",
-% version = "0.3",
-% date = "16 November 2016",
+% version = "0.1.1",
+% date = "1 April 2018",
% codetable = "ISO/ASCII",
% keywords = "LaTeX, goose-article",
-% supported = "http://tdegeus.github.io/GooseLaTeX",
+% supported = "https://github.com/tdegeus/goose-article",
% docstring = "A customized article-class."
% ------------------------------------------------------------------------------
% ==============================================================================
% Basic settings
% ==============================================================================
\NeedsTeXFormat{LaTeX2e}
\ProvidesClass{goose-article}[2016/11/16 v0.3 customized article class]
% load default article class
\LoadClass[a4paper,fleqn,twoside]{article}
% ==============================================================================
% Options
% ==============================================================================
% bibliography [namecite]
% -----------------------
\newif\if@namecite
\let\if@namecite\iffalse
\DeclareOption{namecite}{\let\if@namecite\iftrue}
% Review: [narrow,doublespacing]
% ------------------------------
\newif\if@narrow
\let\if@narrow\iffalse
\DeclareOption{narrow}{\let\if@narrow\iftrue}
\newif\if@doublespacing
\let\if@doublespacing\iffalse
\DeclareOption{doublespacing}{\let\if@doublespacing\iftrue}
+\newif\if@twocolumn
+\let\if@twocolumn\iffalse
+\DeclareOption{twocolumn}{\let\if@twocolumn\iftrue}
+
+% No header
+% ---------
+
+\newif\if@empty
+\let\if@empty\iffalse
+\DeclareOption{empty}{\let\if@empty\iftrue}
+
% Font: [times,garamond,verdana]
% ------------------------------
\newif\if@font
\let\if@font\iffalse
\newif\if@times
\let\if@times\iffalse
\DeclareOption{times}{\let\if@times\iftrue\let\if@font\iftrue}
\newif\if@garamond
\let\if@garamond\iffalse
\DeclareOption{garamond}{\let\if@garamond\iftrue\let\if@font\iftrue}
\newif\if@verdana
\let\if@verdana\iffalse
\DeclareOption{verdana}{\let\if@verdana\iftrue\let\if@font\iftrue}
% Process options
% ---------------
\ProcessOptions\relax
% ==============================================================================
% Packages
% ==============================================================================
\RequirePackage[left=25mm,right=25mm,top=30mm,bottom=30mm]{geometry}
\if@narrow
\newgeometry{left=40mm,right=40mm,top=30mm,bottom=30mm}
\fi
\RequirePackage{fix-cm}
\RequirePackage{graphicx}
\RequirePackage[labelformat=simple]{subcaption}
\RequirePackage{amssymb,amsfonts,bm,cancel}
\RequirePackage{enumerate}
\RequirePackage[fleqn]{amsmath}
\RequirePackage[font={small},labelfont=bf,labelsep=period]{caption}
\RequirePackage[bf]{titlesec}
\RequirePackage[dvipsnames]{xcolor}
\RequirePackage{authblk}
\RequirePackage{hyperref}
+\RequirePackage{multicol}
\renewcommand\thesubfigure{(\alph{subfigure})}
% ==============================================================================
% Citation style
% ==============================================================================
\if@namecite
\RequirePackage{natbib}
\let\oldbibliography\bibliography
\renewcommand{\bibliography}[1]{%
\bibliographystyle{apalike}
\setlength{\bibsep}{3pt plus 0.3ex}
\def\bibfont{\scriptsize}
\oldbibliography{#1}
}
\else
\RequirePackage[square,sort&compress,numbers,comma]{natbib}
\let\oldbibliography\bibliography
\renewcommand{\bibliography}[1]{%
\bibliographystyle{unsrtnat}
\setlength{\bibsep}{3pt plus 0.3ex}
\def\bibfont{\scriptsize}
- \oldbibliography{#1}
+ \begin{multicols}{2}\raggedright\oldbibliography{#1}\end{multicols}
}
\fi
% define DOI command
\newcommand*{\doi}[1]{\sloppy\href{http://dx.doi.org/#1}{doi:~#1}}
-\newcommand*{\eprint}[1]{\sloppy\href{http://arxiv.org/abs/#1}{arXiv:~#1}}
+\newcommand*{\arxivid}[1]{\sloppy\href{http://arxiv.org/abs/#1}{arXiv:~#1}}
% ==============================================================================
% Font
% ==============================================================================
\if@font
\RequirePackage{mathspec}
\RequirePackage{fontspec}
\fi
\if@garamond
\AtBeginDocument{%
\setmathfont(Digits)[Numbers=OldStyle]{Garamond}%
\setromanfont[%
BoldFont={Garamond Bold},%
ItalicFont={Garamond Italic},%
Mapping=tex-text%
]{Garamond}%
}%
\fi
\if@verdana
\AtBeginDocument{%
\defaultfontfeatures{Mapping=tex-text}
\setmainfont[Mapping=tex-text]{Verdana}
\setsansfont{Verdana}
\renewcommand*{\familydefault}{\sfdefault}
\setmathfont(Greek,Digits,Latin){Verdana}
\setmathrm{Verdana}
}%
\fi
\if@times
\AtBeginDocument{%
\setmathfont(Digits)[Numbers=OldStyle]{Times New Roman}%
\setromanfont[%
BoldFont={Times New Roman Bold},%
ItalicFont={Times New Roman Italic},%
Mapping=tex-text%
]{Times New Roman}%
}%
\fi
% ==============================================================================
% Spacing
% ==============================================================================
\if@doublespacing
\AtBeginDocument{%
\usepackage{setspace}
\doublespacing
}%
\fi
+% ==============================================================================
+% Two-columns
+% ==============================================================================
+
+\if@twocolumn
+ \AtBeginDocument{%
+ \twocolumn
+ }%
+\fi
+
% ==============================================================================
% Hyperref
% ==============================================================================
\AtBeginDocument{%
\hypersetup{%
pdftitle=\@title,
citecolor=Blue,
filecolor=black,
linkcolor=Blue,
urlcolor=Blue,
breaklinks,
hidelinks,
bookmarksopen=true,
}%
}
% ==============================================================================
% Headers
% ==============================================================================
% {command}{format}{label}{sep}{before-code}{after-code}
% - section
\titleformat
{\section}
{\normalfont\large\bfseries}
{\thesection}
{1em}
{}
{}
% - subsection
\titleformat
{\subsection}
{\normalfont\normalsize\bfseries}
{\thesubsection}
{1em}
{}
{}
% - paragraph
\titleformat
{\paragraph}
{\normalfont\normalsize\bfseries}
{}
{}
{}
+ {}
% spacing {left}{top}{bottom}, of:
% - section
\titlespacing{\section}
{0pt}{1eM}{0.8em}
% - subsection
\titlespacing{\subsection}
{0pt}{0.8em}{0.4em}
+% - paragraph
+\titlespacing{\paragraph}
+ {0pt}{0.8em}{0em}
% no paragraph indent after header
\AtBeginDocument{%
\makeatletter
\let\orig@afterheading\@afterheading
\def\@afterheading{%
\@afterindentfalse
\orig@afterheading}
\makeatother
}
% ==============================================================================
% Title
% ==============================================================================
% title, authors, affiliation layout
% - contact
\newcommand{\contact}[1]{\date{{#1}}}
% - basic settings
\newcommand{\nl}{\vspace*{-.3em}\protect\\ \fontsize{9}{10.8}\itshape}
\renewcommand\Authfont{\fontsize{11}{14.4}\selectfont}
\renewcommand\Affilfont{\fontsize{9}{10.8}\itshape}
\renewcommand*{\Authsep}{, }
\renewcommand*{\Authand}{, }
\renewcommand*{\Authands}{, }
% - actual layout
\renewcommand{\maketitle}{%
\newpage
\null
\vskip 2em%
\begin{center}%
\let \footnote \thanks
{\Large\bfseries \@title \par}%
\vskip 1.5em%
{\normalsize
\lineskip .5em%
\begin{tabular}[t]{c}%
\@author
\end{tabular}\par}%
\vskip 0.5em%
{\small \@date}%
\end{center}%
\par
\vskip 1.5em
\thispagestyle{fancy}
}
% ==============================================================================
% Header and footer
% ==============================================================================
\RequirePackage{fancyhdr}
\setlength{\headheight}{16pt}
\pagestyle{fancy}
\fancyhead{}
\fancyfoot{}
\renewcommand{\headrulewidth}{0pt}
\renewcommand{\footrulewidth}{0pt}
-\fancyhead[RO,LE]{\scriptsize \thepage}
-\newcommand{\header}[1]{\fancyhead[LO,RE]{\scriptsize\itshape#1}}
+\if@empty
+\else
+\fancyhead[RO,LE]{\thepage}
+\fi
+\newcommand{\header}[1]{\fancyhead[LO,RE]{#1}}
+\newcommand{\headerOdd}[1]{\fancyhead[LO]{#1}}
+\newcommand{\headerEven}[1]{\fancyhead[RE]{#1}}
% ==============================================================================
% Abstract and keywords
% ==============================================================================
\renewenvironment{abstract}{\section*{Abstract}\noindent\ignorespaces}{\ignorespacesafterend}
\newcommand{\keywords}[1]{\vspace*{0.5eM}\noindent\textbf{Keywords: }#1}
diff --git a/docs/theory/readme.pdf b/docs/theory/readme.pdf
index cb45de0..b7802e5 100644
Binary files a/docs/theory/readme.pdf and b/docs/theory/readme.pdf differ
diff --git a/docs/theory/readme.tex b/docs/theory/readme.tex
index 1967279..b9b4735 100644
--- a/docs/theory/readme.tex
+++ b/docs/theory/readme.tex
@@ -1,1171 +1,1192 @@
%!TEX program = xelatex
\documentclass[times,namecite]{goose-article}
\usepackage{lipsum,verbatim,mdframed,array}
\title{%
The Finite Element Method in Continuum Mechanics
}
\author{T.W.J.~de~Geus}
\contact{%
$^*$Contact: %
\href{http://goosefem.geus.me}{goosefem.geus.me}%
\hspace{1mm}--\hspace{1mm} %
\href{mailto:tom@geus.me}{tom@geus.me} %
\hspace{1mm}--\hspace{1mm} %
\href{http://www.geus.me}{www.geus.me}%
}
\hypersetup{pdfauthor={T.W.J. de Geus}}
-\header{%
+\headerEven{\emph{%
The Finite Element Method in Continuum Mechanics, T.W.J.\ de Geus
-}
+}}
\begin{document}
\maketitle
\begin{abstract}
The theory of the Finite Element Method in Continuum Mechanics is discussed in a compact way. This discussion is by no means comprehensive and one is invited to dive in more complete textbooks, but also to get one's hands dirty by implementing some (simple) examples.
\end{abstract}
\setcounter{tocdepth}{2}
\tableofcontents
\section{Statics}
\subsection{The basic concept}
-The discussion starts by considering a static, solid mechanics, problem. Loosely speaking the goal is to find a deformation map, $\vec{x} = \varphi(\vec{X},t)$, that maps a body $\Omega_0$ to a deformed state $\Omega$ that satisfies equilibrium and the boundary conditions applied on $\Gamma$. This is illustrated in Fig.~\ref{fig:problem}, whereby it is emphasized that the body can be subjected to two kinds of boundary conditions:
+The discussion starts by considering a static, solid mechanics, problem. Loosely speaking the goal is to find a deformation map, $\vec{x} = \varphi(\vec{X},t)$, that maps a body $\Omega_0$ to a deformed state $\Omega$ that satisfies equilibrium \emph{and} the boundary conditions applied on $\Gamma$. This is illustrated in Fig.~\ref{fig:problem}, whereby it is emphasised that the body can be subjected to two kinds of boundary conditions:
\begin{itemize}
\item \emph{Essential} or \emph{Dirichlet} boundary conditions on $\Gamma_p$, whereby the displacements are prescribed.
\item \emph{Natural} or \emph{Neumann} boundary conditions on $\Gamma_u$, whereby the tractions are prescribed. Note that traction-free natural boundary conditions are also perfectly acceptable.
\end{itemize}
\begin{figure}[htp]
\centering
\includegraphics[width=.5\textwidth]{figures/problem.pdf}
\caption{Generic problem definition}
\label{fig:problem}
\end{figure}
In practice, one is not looking for the map $\vec{x} = \varphi(\vec{X},t)$, but for the deformation gradient $\bm{F}$, or in fact for a displacement field $\vec{u} = \vec{x} - \vec{X}$. To make things a bit more explicit, the deformation gradient is defined as follows:
\begin{equation}
\vec{x} = \bm{F} \cdot \vec{X}
\end{equation}
hence
\begin{equation}
\bm{F}
=
\frac{\partial \varphi}{\partial \vec{X}}
=
\big( \vec{\nabla}_0 \, \vec{x} \big)^T
=
\bm{I} + \big( \vec{\nabla}_0 \, \vec{u} \big)^T
\end{equation}
\subsection{Momentum balance}
The momentum balance reads
\begin{equation}
\label{eq:static:momentum}
\vec{\nabla} \cdot \bm{\sigma}(\vec{x}) = \vec{0}
\qquad
\vec{x} \in \Omega
\end{equation}
where $\bm{\sigma}$ is the Cauchy stress which depends on the new position $\vec{x}$ and thus on the displacement $\vec{u}$. It has been assumed that all actions are instantaneous (no inertia) and, for simplicity, that there are no body forces. Loosely speaking the interpretation of this equation is that \emph{the sum of all forces vanishes} everywhere in the domain $\Omega$.
The crux of the Finite Element Method is that this non-linear differential equation is solved in a weak sense. I.e.
\begin{equation}
\label{eq:static:int}
\int\limits_\Omega
\vec{\phi}(\vec{x}) \cdot \big[\, \vec{\nabla} \cdot \bm{\sigma}(\vec{x}) \,\big] \;
\mathrm{d}\Omega
=
0
\qquad
\forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d
\end{equation}
where $\vec{\phi}$ are test functions. For reasons that become obvious below, integration by parts is applied, which results in
\begin{equation}
\label{eq:static:weak}
\int\limits_\Omega
\big[\, \vec{\nabla} \vec{\phi}(\vec{x}) \,\big] : \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\int\limits_\Omega
\vec{\nabla}
\cdot
\big[\, \vec{\phi}(\vec{x}) \cdot \bm{\sigma}(\vec{x}) \,\big] \;
\mathrm{d}\Omega
\qquad
\forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d
\end{equation}
In going from Eq.~\eqref{eq:static:int} to Eq.~\eqref{eq:static:weak}, use has been made of the following chain rule
\begin{equation}
\vec{\nabla} \cdot \big[\, \vec{\phi} \cdot \bm{\sigma}^T \,\big]
=
\big[\, \vec{\nabla} \vec{\phi} \,\big] : \bm{\sigma}^T
+
\vec{\phi} \cdot \big[\, \vec{\nabla} \cdot \bm{\sigma} \,\big]
\end{equation}
-together with the symmetry of the Cauchy stress
+Note thereby that the Cauchy stress tensor is symmetric:
\begin{equation}
\bm{\sigma} = \bm{\sigma}^T
\end{equation}
The right-hand side of this equation can be reduced to an area integral by employing Gauss's divergence theorem:
\begin{equation}
- \int\limits_\Omega \vec{\nabla} \cdot \vec{a}(\vec{x}) \; \mathrm{d}\Omega
+ \int\limits_\Omega \vec{\nabla} \cdot \vec{t}(\vec{x}) \; \mathrm{d}\Omega
=
- \int\limits_\Gamma \vec{n}(\vec{x}) \cdot \vec{a}(\vec{x}) \; \mathrm{d}\Gamma
+ \int\limits_\Gamma \vec{n}(\vec{x}) \cdot \vec{t}(\vec{x}) \; \mathrm{d}\Gamma
\end{equation}
where $\vec{n}$ is the normal along the surface $\Gamma$. The result reads
\begin{equation}
\label{eq:static:weak:final}
\int\limits_\Omega
\big[\, \vec{\nabla} \vec{\phi}(\vec{x}) \,\big] : \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\int\limits_\Gamma
\vec{\phi}(\vec{x}) \cdot
\underbrace{
\vec{n}(\vec{x}) \cdot \bm{\sigma}(\vec{x})
}_{
\vec{t}(\vec{x})
} \;
\mathrm{d}\Gamma
\qquad
\forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d
\end{equation}
-\subsection{Discretization}
+\subsection{Discretisation}
-The problem is now discretized using $n$ \emph{nodes} that are connected through \emph{elements}, which define the discretized domain $\Omega^h_0$. Shape functions $N_i(\vec{x})$ are used to extrapolate the nodal quantities throughout the domain $\Omega^h_0$ (and $\Omega^h$). For example the displacement field $\vec{u}(\vec{x})$
+The problem is now discretised using $n$ \emph{nodes} that are connected through \emph{elements}, which define the discretised domain $\Omega^h_0$. Shape functions $N_i(\vec{x})$ are used to interpolate the nodal quantities throughout the domain $\Omega^h_0$ (and $\Omega^h$). For example, the discretisation and interpolation of the displacement field reads
\begin{equation}
\vec{u}(\vec{x},t)
\approx
\vec{u}^h(\vec{x},t)
=
\sum_{m=1}^{n_\mathrm{nodes}} N_m (\vec{x}) \; \vec{u}_m (t)
\end{equation}
whereby we will not further consider the time dependence (denoted with $t$) for now. Following standard Galerkin the test functions are interpolated in the same way, i.e.\
\begin{equation}
-\label{eq:discretization}
+\label{eq:discretisation}
\vec{\phi}(\vec{x})
\approx
\vec{\phi}^h(\vec{x})
=
\sum_{m=1}^{n_\mathrm{nodes}} N_m (\vec{x}) \; \vec{\phi}_m
\end{equation}
-Applied to the problem sketch from Fig.~\ref{fig:problem}, a discretization might look like Fig.~\ref{fig:problem:discretized}. The nodes are clearly marked as circles. The lines connecting the nodes clearly mark the elements which are in this case three-node linear triangles.
+Applied to the problem sketch from Fig.~\ref{fig:problem}, a discretisation might look like Fig.~\ref{fig:problem:discretised}. The nodes are clearly marked as circles. The lines connecting the nodes clearly mark the elements which are in this case three-node linear triangles.
\begin{figure}[htp]
\centering
\includegraphics[width=.25\textwidth]{figures/problem-discretized.pdf}
- \caption{Example of a discretization of the problem in Fig.~\ref{fig:problem} using Finite Elements (linear triangles in this case).}
- \label{fig:problem:discretized}
+ \caption{Example of a discretisation of the problem in Fig.~\ref{fig:problem} using Finite Elements (linear triangles in this case).}
+ \label{fig:problem:discretised}
\end{figure}
Applied to the balance equation the following is obtained
\begin{equation}
\vec{\phi}_m \cdot
\int\limits_{\Omega^h}
\big[\, \vec{\nabla} N_m (\vec{x}) \,\big]
\cdot
\bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\vec{\phi}_m \cdot
\int\limits_{\Gamma^h}
N_m (\vec{x}) \cdot
\vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\qquad
\forall \; \vec{\phi}_m \in \mathbb{R}^d_n
\end{equation}
from which the dependency on $\vec{\phi}_m$ can be dropped:
\begin{equation}
\int\limits_{\Omega^h}
\big[\, \vec{\nabla} N_m (\vec{x}) \,\big]
\cdot
\bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\int\limits_{\Gamma^h}
N_m (\vec{x}) \cdot
\vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\end{equation}
This corresponds to a (non-linear) set of nodal balance equations:
\begin{equation}
\vec{f}_m^\mathrm{int}(\vec{x})
=
\vec{f}_m^\mathrm{ext}(\vec{x})
\end{equation}
with:
\begin{itemize}
\item \emph{Internal forces}
\begin{equation}
\vec{f}_m^\mathrm{int}(\vec{x})
=
\int\limits_{\Omega^h}
\big[\, \vec{\nabla} N_m (\vec{x}) \,\big]
\cdot
\bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
\end{equation}
\item \emph{External forces}
\begin{equation}
\vec{f}_m^\mathrm{ext}(\vec{x})
=
\int\limits_{\Gamma^h}
N_m (\vec{x}) \cdot
\vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\end{equation}
- Note that this term is zero in the interior of the domain, i.e. in $\Omega^h \bigcap \Gamma^h$, while it can be zero or non-zero in $\Gamma^h$ depending on the problem details. Most often it following completely from the boundary conditions of the boundary value problem that we set out to solve.
+ Note that this term is zero in the interior of the domain, i.e. in $\Omega^h \bigcap \Gamma^h$, while it can be zero or non-zero in $\Gamma^h$ depending on the problem details. Most often it follows completely from the boundary conditions of the boundary value problem that we set out to solve.
\end{itemize}
\subsection{Iterative solution -- small strain}
A commonly used strategy to solve the non-linear system is the iterative Newton-Raphson scheme (see Appendix~\ref{sec:newton-raphson}). The idea is thereby to formulate an initial guess for the solution, determine possible residual forces, and use these forces to come to a better guess for the solution. This is continued until the solution has been found, i.e. when the residual vanishes.
-This solution technique is discussed here in the context of small deformations. Assuming the deformations (and therefore rotations) to be small leads to the assumption that $\Omega = \Omega_0$, and thus that $\nabla = \nabla_0$. In this context one typically uses the linear strain tensor
+This solution technique is discussed here in the context of small deformations. Assuming the deformations (and therefore rotations) to be small, leads to the assumption that $\Omega = \Omega_0$, and thus that $\nabla = \nabla_0$. In this context one typically uses the linear strain tensor
\begin{equation}
\bm{\varepsilon}
=
\tfrac{1}{2} \left[ \nabla_0 \vec{u} + \big[\, \nabla_0 \vec{u} \,\big]^T \right]
=
\mathbb{I}_s : \big[\, \nabla_0 \vec{u} \,\big]
\end{equation}
and some (non-linear) relationship between it and the stress as follows
\begin{equation}
\bm{\sigma} = \bm{\sigma} \big( \bm{\varepsilon} \big)
\end{equation}
To further simplify the discussion, the boundary tractions are assumed to be some displacement independent quantity, which is \textit{a-priori} known for all relevant boundary nodes.
The nodal equilibrium equations now reads
\begin{equation}
\label{eq:static:eq-non-lin}
\vec{r}_m(\vec{X}, \vec{u})
=
\vec{f}_m^\mathrm{ext}(\vec{X})
-
\vec{f}_m^\mathrm{int}(\vec{X}, \vec{u})
=
\underline{\vec{0}}
\end{equation}
with
\begin{equation}
\vec{f}_m^\mathrm{int}(\vec{X}, \vec{u})
=
\int\limits_{\Omega^h_0}
\big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big]
\cdot
\bm{\sigma}(\vec{X}, \vec{u}) \;
\mathrm{d}\Omega_0
\end{equation}
The displacement field $\vec{u}$ is now iteratively updated starting from some initial guess, i.e.
\begin{equation}
\vec{u}_{(i+1)} = \vec{u}_{(i)} + \delta \vec{u}
\end{equation}
The update, $\delta \vec{u}$, follows from the linearization of Eq.~\eqref{eq:static:eq-non-lin}. The results of which is
\begin{equation}
\int\limits_{\Omega^h_0}
\big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big]
\cdot
\mathbb{K}\big(\vec{X},\vec{u}_{(i)}\big)
\cdot
\big[\, \vec{\nabla}_0 N_n(\vec{X}) \,\big] \;
\mathrm{d}\Omega_0
\cdot \delta \vec{u}_n
=
\vec{f}_m^\mathrm{ext}(\vec{X})
-
\int\limits_{\Omega^h_0}
\big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big]
\cdot
\bm{\sigma}\big(\vec{X},\vec{u}_{(i)}\big) \;
\mathrm{d}\Omega_0
\end{equation}
where, at each position $\vec{X}$,
\begin{equation}
\mathbb{K}\big(\vec{u}_{(i)}\big)
=
\left. \frac{\partial \bm{\sigma}}{\partial \bm{\varepsilon}} \right|_{\vec{u}_{(i)}}
:
\mathbb{I}_s
\end{equation}
is the \emph{constitutive tangent operator}.
In a shorter notation, the iterative update reads:
\begin{equation}
\bm{K}_{mn,(i)} \cdot \delta \vec{u}_m
=
\vec{f}_m^\mathrm{ext}
-
\vec{f}_{m,(i)}^\mathrm{int}
\end{equation}
with
\begin{equation}
\bm{K}_{mn,(i)}
=
\int\limits_{\Omega^h_0}
\big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big]
\cdot
\mathbb{K}\big(\vec{X}, \vec{u}_{(i)}\big)
\cdot
\big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big] \;
\mathrm{d}\Omega_0
\end{equation}
and
\begin{equation}
\vec{f}_{m,(i)}^\mathrm{int}
=
\int\limits_{\Omega^h_0}
\big[\, \vec{\nabla}_0 N_m(\vec{X}) \,\big]
\cdot
\bm{\sigma}\big(\vec{X}, \vec{u}_{(i)}\big) \;
\mathrm{d}\Omega_0
\end{equation}
\section{Dynamics}
\subsection{Momentum balance}
The next step is to add inertia to the balance of momentum in Eq.~\eqref{eq:static:momentum}:
\begin{equation}
\rho\, \vec{a}
=
\vec{\nabla} \cdot
\bm{\sigma}(\vec{x})
\qquad
\vec{x} \in \Omega
\end{equation}
where $\rho$ is the mass density, and second time derivative of the position $\vec{x}$ is the acceleration $\vec{a} = \ddot{\vec{u}} = \ddot{\vec{x}}$. Note that this function is the continuum equivalent of $\vec{f} = m \vec{a}$.
Like before, this equation is solved in a weak sense
\begin{equation}
\int\limits_\Omega
\rho(\vec{x})\; \vec{\phi}(\vec{x}) \cdot \vec{a} \;
\mathrm{d}\Omega
=
\int\limits_\Omega
\vec{\phi}(\vec{x})
\cdot
\Big[\,
\vec{\nabla}
\cdot
\bm{\sigma}(\vec{x})
\,\Big] \;
\mathrm{d}\Omega
\qquad
\forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d
\end{equation}
Integration by parts results in
\begin{equation}
\int\limits_\Omega
\rho(\vec{x})\; \vec{\phi}(\vec{x}) \cdot \vec{a} \;
\mathrm{d}\Omega
=
\int\limits_\Gamma
\vec{\phi}(\vec{x}) \cdot \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
-
\int\limits_\Omega
\big[\, \vec{\nabla} \vec{\phi}(\vec{x}) \,\big]
:
\bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
\qquad
\forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d
\end{equation}
-which is discretized as before:
+which is discretised as before:
\begin{equation}
\vec{\phi}_m \cdot
\int\limits_\Omega
\rho(\vec{x})\; N_m(\vec{x})\; N_n(\vec{x}) \;
\mathrm{d}\Omega \;
\vec{a}_n
=
\vec{\phi}_m \cdot
\int\limits_\Gamma
N_m(\vec{x})\; \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
-
\vec{\phi}_m \cdot
\int\limits_\Omega
\big[\, \vec{\nabla} N_m(\vec{x}) \,\big]
:
\bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
- \qquad
- \forall \; \vec{\phi}_m \in \mathbb{R}^d_n
\end{equation}
-This is independent of the test functions, hence:
+which should still hold for all $\vec{\phi}_m \in \mathbb{R}^d_n$. However, this dependence can now trivially be dropped. Hence:
\begin{equation}
\label{eq:dynamics:system}
\underbrace{
\int\limits_\Omega
\rho(\vec{x})\; N_m(\vec{x})\; N_n(\vec{x}) \;
\mathrm{d}\Omega
}_{\displaystyle
M_{mn}(\vec{x})
} \;
\vec{a}_n
=
\underbrace{
\int\limits_\Gamma
N_m(\vec{x})\; \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
}_{\displaystyle
\vec{f}_m^\mathrm{ext}(\vec{x})
}
-
\underbrace{
\int\limits_\Omega
\big[\, \vec{\nabla} N_m(\vec{x}) \,\big]
:
\bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
}_{\displaystyle
\vec{f}_m^\mathrm{int}(\vec{x})
}
\end{equation}
This is denoted as follows
\begin{equation}
M_{mn}(\vec{x})\; \vec{a}_n
=
\vec{f}_m^\mathrm{ext}(\vec{x})
-
\vec{f}_m^\mathrm{int}(\vec{x})
\end{equation}
Where $M_{mn}(\vec{x})$ is the \emph{mass matrix}, $\vec{f}_m^\mathrm{ext}(\vec{x})$ are the \emph{external forces}, and $\vec{f}_m^\mathrm{int}(\vec{x})$ are the \emph{internal forces}.
-\subsection{Time discretization}
+\subsection{Time discretisation}
\label{eq:dynamics:time}
-To solve the second order differential equation in time, one typically also discretizes time with some finite difference based scheme. For this we also need to distinguish the velocity field $\vec{v} = \dot{\vec{u}} = \dot{\vec{x}}$.
+To solve the second order differential equation in time, one typically also discretises time with some finite difference based scheme. For this we also need to distinguish the velocity field $\vec{v} = \dot{\vec{u}} = \dot{\vec{x}}$.
\subsubsection{Verlet}
-In the case that the forces are velocity independent one can use the Velocity Verlet scheme, in which energy is conserved, i.e.\ it is reversible in time. The scheme reads are follows:
+In the case that the forces are velocity independent one can use the Velocity Verlet scheme, in which energy is conserved, i.e.\ it is reversible in time. The scheme reads as follows:
\begin{enumerate}
\item Compute the velocity on a dummy time grid:
\begin{equation}
\vec{v}_{n+1/2} = \vec{v}_{n-1/2} + \Delta_t \; \vec{a}_n
\end{equation}
Note that this involves solving Eq.~\eqref{eq:dynamics:system} for $\vec{a}_n$. From this it is directly obvious why the forces need to be velocity independent.
\item Update the positions
\begin{equation}
\vec{u}_{n+1} = \vec{u}_n + \Delta_t \;\vec{a} \vec{v}_{n + 1/2}
\end{equation}
\end{enumerate}
\subsubsection{Velocity Verlet}
In the case that the forces are dependent on the velocity (i.e.\ when there is damping) the previous scheme cannot be used anymore. An additional approximation has to be made in this case. Note that this is not very problematic as energy is in no case conserved. The resulting scheme reads:
\begin{enumerate}
\item Compute the position at $t_{n+1} = t_{n} + \Delta_t$:
\begin{equation}
\vec{u}_{n+1}
=
\vec{u}_{n} + \Delta_t \vec{v}_{n} + \tfrac{1}{2} \Delta_t^2 \vec{a}_{n}
\end{equation}
\item Estimate the velocity at $t_{n+1} = t_{n} + \Delta_t$ (involves solving Eq.~\eqref{eq:dynamics:system}):
\begin{equation}
\hat{\vec{v}}_{n+1}
=
\vec{v}_{n}
+
\tfrac{1}{2} \Delta_t \Big[\,
\vec{a}_{n} + \vec{a} ( \vec{u}_{n+1} , \vec{v}_{n} + \Delta_t \vec{a}_{n} , t_{n+1} ) \,
\Big]
\end{equation}
\item Correct $\hat{\vec{v}}_{n+1}$ (involves solving Eq.~\eqref{eq:dynamics:system}):
\begin{equation}
\vec{v}_{n+1}
=
\vec{v}_{n}
+
\tfrac{1}{2} \Delta_t \Big[\,
\vec{a}_{n} + \vec{a} ( \vec{u}_{n+1} , \hat{\vec{v}}_{n+1} , t_{n+1} ) \,
\Big]
\end{equation}
\end{enumerate}
\subsection{Approximations}
For problems where the local volume is conversed (either weakly through slightly compressible elasticity, or strongly by adding an incompressibility constraint) it makes sense to assume the mass matrix constant, as any change of volume results in an equivalent change of the density. In that case
\begin{equation}
\int\limits_{\Omega}
\rho(\vec{x})
\;\mathrm{d}\Omega
=
\int\limits_{\Omega_0}
\rho(\vec{X})
\;\mathrm{d}\Omega_0
\end{equation}
Which results in:
\begin{equation}
M_{mn}(\vec{x})
=
\int\limits_{\Omega_0}
\rho(\vec{X})\; N_m(\vec{X})\; N_n(\vec{X}) \;
\mathrm{d}\Omega_0
=
\mathrm{constant}
\end{equation}
To enhance computational efficiency, it may be a good option concentrate the mass to point masses on the nodes. This has to strong advantage that the mass matrix becomes diagonal, known as the \emph{lumped mass matrix}. Consequently, instead of solving a linear system one just has to solve fully independent equations.
\appendix
\section{Nomenclature}
\begin{itemize}
\item Vector
\begin{equation}
\vec{u} = \sum\limits_i u_i \vec{e}_i
\end{equation}
for example for Cartesian coordinates
\begin{equation}
\vec{u} = u_x \vec{e}_x + u_y \vec{e}_y + u_z \vec{e}_z
\end{equation}
%
\item Second-order tensor
\begin{equation}
\bm{A} = \sum\limits_i \sum\limits_j A_{ij} \vec{e}_i \vec{e}_j
\end{equation}
%
\item Fourth-order tensor
\begin{equation}
\mathbb{A} = \sum\limits_i \sum\limits_j \sum\limits_k \sum\limits_l A_{ijkl} \vec{e}_i \vec{e}_j \vec{e}_k \vec{e}_l
\end{equation}
%
\item Divergence
\begin{equation}
\vec{\nabla} \cdot \bm{\sigma} = \frac{ \partial \sigma_{ij} }{ \partial x_i }
\end{equation}
%
\item Double tensor contraction
\begin{equation}
C = \bm{A} : \bm{B} = A_{ij} B_{ji}
\end{equation}
\end{itemize}
\section{Shape functions}
-In the Finite Element Method a geometry is discretized using nodes. The nodes are grouped in elements which define the domain $\Omega^h$. The crux of the method is that nodal quantities, for example $\vec{u}_i$, are interpolated throughout the discretized domain $\Omega^h$ using shape functions $N_i (\vec{x})$. The shape functions are polynomials that spanned by the nodes, such that they constitute a partition of unity and satisfy $N_i (\vec{x}_j) = \delta_{ij}$, i.e. it is one in the node $i$, and zero in all other nodes. This implies that, while each shape function is globally supported, $N_i (\vec{x}) \neq 0$ only in the elements containing node $i$. For all other elements it is zero (even in between the nodes). This again implies that $\partial N_i / \partial \vec{x}$ is discontinuous across element boundaries.
+In the Finite Element Method a geometry is discretised using nodes. The nodes are grouped in elements which define the domain $\Omega^h$. The crux of the method is that nodal quantities, for example $\vec{u}_m$, are interpolated throughout the discretised domain $\Omega^h$ using shape functions $N_m (\vec{x})$. The shape functions are polynomials that spanned by the nodes, such that they constitute a partition of unity and satisfy $N_m (\vec{x}_n) = \delta_{mn}$, i.e. it is one in the node $i$, and zero in all other nodes. This implies that, while each shape function is globally supported, $N_m (\vec{x}) \neq 0$ only in the elements containing node $i$. For all other elements it is zero (even in between the nodes). This again implies that $\partial N_m / \partial \vec{x}$ is discontinuous across element boundaries.
-For a one-dimensional problem comprising four linear elements and five nodes the shape functions are sketched in Fig.~\ref{fig:shape:1d-domain}. To emphasize the global support of each shape function the shape function of node 2 is isolated in Fig.~\ref{fig:shape:1d-domain-2}. As can be seen, node 2 is only non-zero in elements 1 and 2, while it is zero everywhere else. To evaluate for example $f_2$, the internal force one node 2, the integration can be limited to these two elements:
+For a one-dimensional problem comprising four linear elements and five nodes, the shape functions are sketched in Fig.~\ref{fig:shape:1d-domain}. To emphasise the global support of each shape function the shape function of node 2 is isolated in Fig.~\ref{fig:shape:1d-domain-2}. As can be seen, node 2 is only non-zero in elements 1 and 2, while it is zero everywhere else. To evaluate for example $f_2$, the internal force one node 2, the integration can be limited to these two elements:
\begin{equation}
f_2
=
\int\limits_{\Omega^{(1)}}
\frac{\mathrm{d} N^{(1)}_2}{\mathrm{d} x} \;
\sigma (x) \;
\mathrm{d}\Omega^{(1)}
+
\int\limits_{\Omega^{(2)}}
\frac{\mathrm{d} N^{(2)}_2}{\mathrm{d} x} \;
\sigma (x) \;
\mathrm{d}\Omega^{(2)}
\end{equation}
\begin{figure}[htp]
\centering
\captionsetup[subfigure]{justification=centering}
\begin{minipage}[t]{.5\textwidth}
\centering
\includegraphics[width=1.\textwidth]{figures/shape-functions-1d.pdf}
- \subcaption{All shape functions $N_i(x)$, each with a different color.}
+ \subcaption{All shape functions $N_m(x)$, each with a different colour.}
\label{fig:shape:1d-domain}
\end{minipage}
\\
\begin{minipage}[t]{.5\textwidth}
\centering
\includegraphics[width=1.\textwidth]{figures/shape-functions-1d-node-2.pdf}
\subcaption{Isolated shape function $N_2(x)$.}
\label{fig:shape:1d-domain-2}
\end{minipage}
- \caption{(a) Shape functions for a one-dimensional domain discretized using four two node, linear, element. (b) Isolated shape functions for node 2. The node numbers are displayed using a unique color, while the element numbers are shown in black in a different font, in between the nodes.}
+ \caption{(a) Shape functions for a one-dimensional domain discretised using four two node, linear, element. (b) Isolated shape functions for node 2. The node numbers are displayed using a unique colour, while the element numbers are shown in black in a different font, in between the nodes.}
\label{fig:shape:1d-domain:fig}
\end{figure}
-By now it should be clear that the above allows us assemble $\underline{f}$ element-by-element. For the example in Fig.~\ref{fig:shape:1d-domain:fig} this is graphically shown in Fig.~\ref{fig:shape:1d-domain:assembly}.
+By now it should be clear that the above allows the global system $f_m$ to be assembled element-by-element. For the example in Fig.~\ref{fig:shape:1d-domain:fig} this is graphically shown in Fig.~\ref{fig:shape:1d-domain:assembly} where the indices show that the \emph{shape functions} are evaluated compared to some generic element definition.
\begin{figure}[htp]
\centering
\captionsetup[subfigure]{justification=centering}
\begin{minipage}[t]{.5\textwidth}
\centering
\includegraphics[width=1.\textwidth]{figures/shape-functions-1d-element-0.pdf}
\end{minipage}
\\
$+$
\\
\begin{minipage}[t]{.5\textwidth}
\centering
\includegraphics[width=1.\textwidth]{figures/shape-functions-1d-element-1.pdf}
\end{minipage}
\\
$+$
\\
\begin{minipage}[t]{.5\textwidth}
\centering
\includegraphics[width=1.\textwidth]{figures/shape-functions-1d-element-2.pdf}
\end{minipage}
\\
$+$
\\
\begin{minipage}[t]{.5\textwidth}
\centering
\includegraphics[width=1.\textwidth]{figures/shape-functions-1d-element-3.pdf}
\end{minipage}
\caption{Element-by-element assembly of the problem in Fig.~\ref{fig:shape:1d-domain:fig}.}
\label{fig:shape:1d-domain:assembly}
\end{figure}
-
-where the indices show that the \emph{shape functions} are evaluated compared to some generic element definition.
-
\section{Iso-parametric transformation and quadrature}
A very important concept in the Finite Element Method is the iso-parametric transformation. It allows an arbitrarily shaped element with volume $\Omega^e$ to be mapped onto a generic \emph{iso-parametric element} of constant volume $Q$, see Fig.~\ref{fig:isoparametric}. Using this mapping it is easy to perform interpolation and numerical quadrature using a generic implementation.
\begin{figure}[htp]
\centering
\includegraphics[width=.5\textwidth]{figures/isoparametric-transform.pdf}
- \caption{Iso-parameteric transformation for a four node, bi-linear, quadrilateral element.}
+ \caption{Iso-parametric transformation for a four node, bi-linear, quadrilateral element.}
\label{fig:isoparametric}
\end{figure}
The mapping between the generic domain $Q$ and the physical domain $\Omega^e$ is as follows
\begin{equation}
\label{eq:isoparametric:coordinate}
- \vec{x} ( \vec{\xi} ) = N_m^{e} \; \underline{x}_m^e
+ \vec{x} ( \vec{\xi} ) = N_m^{e} \; \vec{x}_m^e
\end{equation}
-where the column $\underline{x}^e$ contains the position vectors of the element nodes in the actual coordinates (i.e.\ spanning $\Omega^e$). In order to perform the quadrature on $Q$ we must also map the gradient operator:
+where the column $\vec{x}_m^e$ contains the position vectors of the element nodes in the actual coordinates (i.e.\ spanning $\Omega^e$). In order to perform the quadrature on $Q$ we must also map the gradient operator:
\begin{equation}
\vec{\nabla}_{\xi}\,
=
\vec{e}_i \frac{\partial}{\partial \xi_i}
=
\vec{e}_i \frac{\partial x_j(\vec{\xi})}{\partial \xi_i} \frac{\partial}{\partial x_j}
=
\vec{e}_i \frac{\partial x_j(\vec{\xi})}{\partial \xi_i} \vec{e}_j \cdot \vec{e}_k \frac{\partial}{\partial x_k}
=
\big[\, \vec{\nabla}_{\xi}\, \vec{x}(\vec{\xi}) \,\big] \cdot \vec{\nabla}
=
\bm{J}(\vec{\xi}) \cdot \vec{\nabla}
\end{equation}
i.e.
\begin{equation}
J_{ij} = \frac{\partial x_i}{\partial \xi_j}
\end{equation}
from which the inverse relationship trivially follows
\begin{equation}
\vec{\nabla} = \bm{J}^{-1}(\vec{\xi}) \cdot \vec{\nabla}_{\xi}\,
\end{equation}
where $\bm{J}$ is the Jacobian which can be obtained from the gradient of the shape functions with respect to the iso-parametric coordinates and the nodal coordinates in the actual coordinates through Eq.~\eqref{eq:isoparametric:coordinate}:
\begin{equation}
\bm{J}(\vec{\xi})
=
\vec{\nabla}_{\xi}\, \vec{x}(\vec{\xi})
=
\vec{\nabla}_{\xi}\, \big[\, N_m^{e} \; \vec{x}_m^e \,\big]
=
\big[\, \vec{\nabla}_{\xi}\, N_m^{e} \,\big] \; \vec{x}_m^e
\end{equation}
The gradient of the shape functions with respect to the actual coordinates can now be computed though
\begin{equation}
\vec{\nabla} N_m^{e}
=
\bm{J}^{-1}(\vec{\xi}) \cdot \big[\, \vec{\nabla}_{\xi}\, N_m^{e} \,\big]
\end{equation}
The change of variables should also be accounted for in the integration:
\begin{equation}
\int_\Omega (\ldots) \; \mathrm{d} \Omega
=
\iiint (\ldots) \; \mathrm{d} x \, \mathrm{d} y \, \mathrm{d} z
=
\iiint (\ldots) \; \left| \frac{\partial \left( x, y, z \right)}{\partial \left( \xi, \eta, \mu \right)} \right| \mathrm{d} \xi \, \mathrm{d} \eta \, \mathrm{d} \mu
\end{equation}
where the determinant of the Jacobian
\begin{equation}
\frac{\partial \left( x, y, z \right)}{\partial \left( \xi, \eta, \mu \right)}
=
\left|
\begin{array}{ccc}
\displaystyle \frac{\partial x}{\partial \xi } &
\displaystyle \frac{\partial x}{\partial \eta} &
\displaystyle \frac{\partial x}{\partial \mu }
\\
\displaystyle \frac{\partial y}{\partial \xi } &
\displaystyle \frac{\partial y}{\partial \eta} &
\displaystyle \frac{\partial y}{\partial \mu }
\\
\displaystyle \frac{\partial z}{\partial \xi } &
\displaystyle \frac{\partial z}{\partial \eta} &
\displaystyle \frac{\partial z}{\partial \mu }
\end{array}
\right|
\end{equation}
We can now make use of the Jacobian that we already obtained
\begin{equation}
\int_\Omega (\ldots) \; \mathrm{d} \Omega
=
\iiint (\ldots) \det \left( \bm{J}(\vec{\xi}) \right) \; \mathrm{d} \xi \, \mathrm{d} \eta \, \mathrm{d} \mu
=
\int_Q (\ldots) \det \left( \bm{J}(\vec{\xi}) \right) \; \mathrm{d} Q
\end{equation}
Note that another way to obtain this result is using
\begin{align}
\mathrm{d} \Omega
&=
\mathrm{d} \vec{x}_0 \times \mathrm{d} \vec{x}_1 \cdot \mathrm{d} \vec{x}_2
=
\left[ \mathrm{d} \vec{x}_0 \cdot \bm{J}(\vec{\xi}) \right] \times
\left[ \mathrm{d} \vec{x}_1 \cdot \bm{J}(\vec{\xi}) \right] \cdot
\left[ \mathrm{d} \vec{x}_2 \cdot \bm{J}(\vec{\xi}) \right]
=
\det \big( \bm{J}(\vec{\xi}) \big)\,
\mathrm{d} \vec{\xi}_0 \times \mathrm{d} \vec{\xi}_1 \cdot \mathrm{d} \vec{\xi}_2
\\
&=
\det \big( \bm{J}(\vec{\xi}) \big)\, \mathrm{d} Q
\end{align}
-Numerical quadrature can be formulated (exactly) on the master element. It corresponds to taking the weighted sum of the integrand evaluated at specific \emph{quadrature points} (or \emph{integration points}).
-
-The practical implications, for example for the internal force, are as follows:
+Numerical quadrature can be formulated (exactly) on the master element. It corresponds to taking the weighted sum of the integrand evaluated at specific \emph{quadrature points} (or \emph{integration points}). The practical implications, for example for the internal force, are as follows:
\begin{align}
\vec{f}_m^e
&=
\int\limits_{\Omega^e}
\big[\, \vec{\nabla} N_m \,\big]
\cdot
\bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\int\limits_{Q}
\big[\, \vec{\nabla} N_m \,\big]
\cdot
\bm{\sigma}(\vec{x}) \;
\det \big( \bm{J}(\vec{\xi}) \big) \;
\mathrm{d}Q
\\
&=
\sum_{k}^{n_k}
w_k
\big[\, \vec{\nabla} N_m \,\big]_{\vec{\xi} = \vec{\xi}_k}
\cdot
\bm{\sigma}\big(\vec{x}(\vec{\xi}_k)\big) \;
\det \big( \bm{J}(\vec{\xi}_k) \big) \;
\end{align}
\subsection*{Total Lagrange}
To obtain
\begin{equation}
- \vec{x}(\vec{\xi}),\, \vec{\nabla}_0,\,\, \text{and}\, \int\limits_{\Omega_0} . \;\mathrm{d}\Omega
+ \vec{x}(\vec{\xi}),\, \vec{\nabla}_0,\,\, \text{and}\, \int\limits_{\Omega_0} (\ldots) \;\mathrm{d}\Omega
\end{equation}
simply replace $\vec{x}_m^e$ with $\vec{X}_m^e$ in Eq.~\eqref{eq:isoparametric:coordinate}. For this reason the same element implementation can be used in small strain and finite strain (total Lagrange and updated Lagrange), by providing either $\vec{X}_m^e$ or $\vec{X}_m^e + \vec{u}_m^e$ as input.
\section{Newton-Raphson in one dimension}
\label{sec:newton-raphson}
The objective is to find $x$ such that
\begin{equation}
r(x) = 0
\end{equation}
An initial guess is made for $x$ which is (hopefully) iteratively improved. This iterative value is denoted using $x_{(i)}$. It is updated by making use of the following Taylor expansion
\begin{equation}
r \big( x_{(i+1)} \big)
=
r \big( x_{(i)} \big)
+
\left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \delta x
+
\mathcal{O} \big( \delta x^2 \big)
\approx
0
\end{equation}
where
\begin{equation}
\delta x = x_{(i+1)} - x_{(i)}
\end{equation}
The value of $\delta x$ can be obtained by neglecting higher order terms:
\begin{equation}
r \big( x_{(i)} \big)
+
\left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \delta x
=
0
\end{equation}
From which it follows that
\begin{equation}
\delta x
=
- \left[ \left. \frac{\mathrm{d} r}{\mathrm{d} x} \right|_{x = x_{(i)}} \right]^{-1}
r \big( x_{(i)} \big)
\end{equation}
Thereafter, the improved `guess' of the solution is
\begin{equation}
x_{(i+1)} = x_{(i)} + \delta x
\end{equation}
These steps are then repeated, until the solution has been reached within a certain accuracy $\epsilon$:
\begin{equation}
\left| r \big( x_{(i+1)} \big) \right| < \epsilon
\end{equation}
-
The iterative scheme is well understood from the illustration in Fig.~\ref{fig:newton-raphson}.
\begin{figure}[htp]
\centering
- \includegraphics[width=.3\textwidth]{figures/newton-raphson.pdf}
+ \includegraphics[width=.4\textwidth]{figures/newton-raphson.pdf}
\caption{Schematic of the Newton-Raphson iterations in a one-dimensional problem.}
\label{fig:newton-raphson}
\end{figure}
\section{B-matrix}
\subsection{Discretisation}
-The discretisation from Eq.~\eqref{eq:discretization} can be written more precisely as
+The discretisation from Eq.~\eqref{eq:discretisation} can be written more precisely as
\begin{equation}
\phi_x(\vec{x}) \vec{e}_x +
\phi_y(\vec{x}) \vec{e}_y +
\phi_z(\vec{x}) \vec{e}_z
\approx
\sum_{m=1}^{n_\mathrm{nodes}}
N_m (\vec{x}) \; \phi_{m,x} \vec{e}_x +
N_m (\vec{x}) \; \phi_{m,y} \vec{e}_y +
N_m (\vec{x}) \; \phi_{m,z} \vec{e}_z
\end{equation}
whereby $\vec{\phi}(\vec{x})$ is a vector field such as the test functions, the displacement, the forces, etc. The resulting shape function gradients then read
\begin{align}
&
\vec{\nabla} \vec{\phi}(\vec{x})
\approx
\sum_{m=1}^{n_\mathrm{nodes}}
\bigg(
\vec{e}_x \frac{\partial}{\partial x} +
\vec{e}_y \frac{\partial}{\partial y} +
\vec{e}_z \frac{\partial}{\partial z}
\bigg)
\bigg(
N_m (\vec{x}) \; \phi_{m,x} \vec{e}_x +
N_m (\vec{x}) \; \phi_{m,y} \vec{e}_y +
N_m (\vec{x}) \; \phi_{m,z} \vec{e}_z
\bigg)
\\
&
=
\phi_{m,x} \frac{\partial N_m}{\partial x} \vec{e}_x \vec{e}_x + \nonumber
\phi_{m,y} \frac{\partial N_m}{\partial x} \vec{e}_x \vec{e}_y + \nonumber
\phi_{m,z} \frac{\partial N_m}{\partial x} \vec{e}_x \vec{e}_z
\\
&
+
\phi_{m,x} \frac{\partial N_m}{\partial y} \vec{e}_y \vec{e}_x +
\phi_{m,y} \frac{\partial N_m}{\partial y} \vec{e}_y \vec{e}_y +
\phi_{m,z} \frac{\partial N_m}{\partial y} \vec{e}_y \vec{e}_z
\\
&
+
\phi_{m,x} \frac{\partial N_m}{\partial z} \vec{e}_z \vec{e}_x +
\phi_{m,y} \frac{\partial N_m}{\partial z} \vec{e}_z \vec{e}_y +
\phi_{m,z} \frac{\partial N_m}{\partial z} \vec{e}_z \vec{e}_z
\end{align}
Which can denote in short as
\begin{equation}
\label{eq:bmatrix:discretisation}
\vec{\nabla} \vec{\phi} (\vec{x})
\approx
\sum_{m=1}^{n_\mathrm{nodes}}
\mathcal{B}_m (\vec{x}) \cdot \vec{\phi}_m
\end{equation}
whereby $\mathcal{B}_m$ is third order tensor that collects the gradient of the shape function of a node $m$. The only non-zero components are
\begin{alignat}{3}
\mathcal{B}_{m,x x x} &= \frac{\partial N_m}{\partial x} & \qquad
\mathcal{B}_{m,y x x} &= \frac{\partial N_m}{\partial y} & \qquad
\mathcal{B}_{m,z x x} &= \frac{\partial N_m}{\partial z}
\nonumber
\\
\mathcal{B}_{m,x y y} &= \frac{\partial N_m}{\partial x} & \qquad
\mathcal{B}_{m,y y y} &= \frac{\partial N_m}{\partial y} & \qquad
\mathcal{B}_{m,z y y} &= \frac{\partial N_m}{\partial z}
\\
\mathcal{B}_{m,x z z} &= \frac{\partial N_m}{\partial x} & \qquad
\mathcal{B}_{m,y z z} &= \frac{\partial N_m}{\partial y} & \qquad
\mathcal{B}_{m,z z z} &= \frac{\partial N_m}{\partial z}
\end{alignat}
\subsection{Applied to the weak form}
By means of example the organisation using the `B-matrix' can be applied to the weak form of Eq.~\eqref{eq:static:weak:final}, which reads
\begin{equation}
\int\limits_\Omega
\big[\, \vec{\nabla} \vec{\phi}(\vec{x}) \,\big] : \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\int\limits_\Gamma
\vec{\phi}(\vec{x}) \cdot
\vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\qquad
\forall \; \vec{\phi}(\vec{x}) \in \mathbb{R}^d
\end{equation}
Substituting Eq.~\eqref{eq:bmatrix:discretisation} results in
\begin{equation}
\int\limits_{\Omega^h}
\left( \mathcal{B}_m \cdot \vec{\phi}_m \right) : \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\int\limits_{\Gamma^h}
N_m (\vec{x}) \; \vec{\phi}_m \cdot \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\qquad
\forall \; \vec{\phi}_m \in \mathbb{R}^d_n
\end{equation}
from which we can take the nodal values of the test functions $\vec{\phi}_m$, that no longer depend on space, out of the volume integral:
\begin{equation}
\vec{\phi}_m \cdot
\int\limits_{\Omega^h}
\mathcal{B}_m^T : \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\vec{\phi}_m \cdot
\int\limits_{\Gamma^h}
N_m(\vec{x}) \; \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\qquad
\forall \; \vec{\phi}_m \in \mathbb{R}^d_n
\end{equation}
From which it is obvious that the dependence on $\vec{\phi}_m$ can be dropped, resulting in the discretised weak form:
\begin{equation}
\label{eq:bmatrix:weak_form:discrete}
\int\limits_{\Omega^h}
\mathcal{B}_m^T : \bm{\sigma}(\vec{x}) \;
\mathrm{d}\Omega
=
\int\limits_{\Gamma^h}
N_m(\vec{x}) \; \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\end{equation}
\subsection{Applied to linear elasticity}
As a further illustration linear elasticity can be assumed. Application to Eq.~\eqref{eq:bmatrix:weak_form:discrete} trivially results in
\begin{equation}
\int\limits_{\Omega_0^h}
\tilde{\mathcal{B}}_m^T : \Big( \mathbb{C} : \bm{\varepsilon}(\vec{x}) \Big) \;
\mathrm{d}\Omega
=
\int\limits_{\Gamma^h}
N_m(\vec{x}) \; \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\end{equation}
where $\tilde{\cdot}$ emphasises that $\tilde{\mathcal{B}}_m$ are the gradients of the shape functions with respect to the coordinates in the, undeformed, reference frame. Now, because $\mathbb{C}$ symmetrises, this expression can be written in terms of the displacement gradient directly:
\begin{equation}
\int\limits_{\Omega_0^h}
\tilde{\mathcal{B}}_m^T : \Big( \mathbb{C} : \vec{\nabla}\vec{u}(\vec{x}) \Big) \;
\mathrm{d}\Omega
=
\int\limits_{\Gamma^h}
N_m(\vec{x}) \; \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\end{equation}
The displacement gradient is evaluated discretely using Eq.~\eqref{eq:bmatrix:weak_form:discrete}, resulting in:
\begin{equation}
\int\limits_{\Omega_0^h}
\tilde{\mathcal{B}}_m^T : \Big( \mathbb{C} : ( \tilde{\mathcal{B}}_n \cdot \vec{u}_n) \Big) \;
\mathrm{d}\Omega
=
\int\limits_{\Gamma^h}
N_m(\vec{x}) \; \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\end{equation}
whereby, again, the nodal displacements, $\vec{u}_n$, are not a function of the position and can be taken out of the integral, such as to finally have a system that can be solved:
\begin{equation}
\int\limits_{\Omega_0^h}
\tilde{\mathcal{B}}_m^T : \Big( \mathbb{C} : \tilde{\mathcal{B}}_n \Big) \;
\mathrm{d}\Omega
\cdot \vec{u}_n
=
\int\limits_{\Gamma^h}
N_m(\vec{x}) \; \vec{t}(\vec{x}) \;
\mathrm{d}\Gamma
\end{equation}
\section{Cylindrical coordinates}
-\begin{figure}[htp]
- \centering
- \includegraphics[width=.3\textwidth]{figures/cylindrical_coordinates}
- \caption{Cylindrical coordinate system spanned by $(\vec{e}_r, \vec{e}_\theta, \vec{e}_z)$.}
- \label{fig:cylindrical}
-\end{figure}
+\paragraph{Coordinates}
The cylindrical coordinate system is spanned by $(\vec{e}_r, \vec{e}_\theta, \vec{e}_z)$ as visualised in Fig.~\ref{fig:cylindrical}. From this figure, the transition from cylindrical coordinates to Cartesian coordinates trivially follows to read
\begin{equation}
x = r \cos \theta \qquad
y = r \sin \theta \qquad
z = z
\end{equation}
and the inverse relation to read
\begin{equation}
r = \sqrt{ x^2 + y^2 } \qquad
\theta = \mathrm{arctan}\left( y/x \right) \qquad
z = z
\end{equation}
Likewise, the unit vectors $(\vec{e}_r, \vec{e}_\theta, \vec{e}_z)$ can be expressed in terms of the unit vectors that span the Cartesian basis:
\begin{align}
\vec{e}_r &= \cos \theta \; \vec{e}_x + \sin \theta \; \vec{e}_y \nonumber \\
\vec{e}_\theta &= - \sin \theta \; \vec{e}_x + \cos \theta \; \vec{e}_y \nonumber \\
\vec{e}_z &= \vec{e}_z \label{eq:cylindrical:unit}
\end{align}
-
The gradient operator for cylindrical coordinates reads
\begin{equation}
\vec{\nabla} =
\vec{e}_r \frac{\partial}{\partial r} +
\frac{1}{r} \vec{e}_\theta \frac{\partial}{\partial \theta} +
\vec{e}_z \frac{\partial}{\partial z}
\end{equation}
whereby from Eq.~\eqref{eq:cylindrical:unit} the derivates of the unit vector follow to be
\begin{alignat}{3}
\frac{\partial \vec{e}_r}{\partial r } &= 0 &\qquad \frac{\partial \vec{e}_\theta}{\partial r } &= 0 &\qquad \frac{\partial \vec{e}_z}{\partial r } &= 0 \nonumber \\
\frac{\partial \vec{e}_r}{\partial \theta} &= \vec{e}_\theta &\qquad \frac{\partial \vec{e}_\theta}{\partial \theta} &= - \vec{e}_r &\qquad \frac{\partial \vec{e}_z}{\partial \theta} &= 0 \nonumber \\
\frac{\partial \vec{e}_r}{\partial z } &= 0 &\qquad \frac{\partial \vec{e}_\theta}{\partial z } &= 0 &\qquad \frac{\partial \vec{e}_z}{\partial z } &= 0
\end{alignat}
+\begin{figure}[htp]
+ \centering
+ \includegraphics[width=.3\textwidth]{figures/cylindrical_coordinates}
+ \caption{Cylindrical coordinate system spanned by $(\vec{e}_r, \vec{e}_\theta, \vec{e}_z)$.}
+ \label{fig:cylindrical}
+\end{figure}
+
+\paragraph{Discretisation}
+
A discretised vector field $\vec{\phi}(\vec{x})$ reads
\begin{equation}
\phi_r (\vec{x}) \vec{e}_r + \phi_\theta (\vec{x}) \vec{e}_\theta + \phi_z (\vec{x}) \vec{e}_z
\approx
\sum\limits_{m=1}^{n_\mathrm{nodes}}
N^m (\vec{x}) \phi^{m}_{r} \vec{e}_r +
N^m (\vec{x}) \phi^{m}_{\theta} \vec{e}_\theta +
N^m (\vec{x}) \phi^{m}_{z} \vec{e}_z
\end{equation}
And consequently its discretised gradient reads
\begin{align}
&\vec{\nabla} \vec{\phi} (\vec{x})
\approx
\sum\limits_{m=1}^{n_\mathrm{nodes}}
\Bigg[
\bigg(
\vec{e}_r \frac{\partial}{\partial r} +
\frac{1}{r} \vec{e}_\theta \frac{\partial}{\partial \theta} +
\vec{e}_z \frac{\partial}{\partial z}
\bigg)
\bigg(
- N^m (\vec{x}) \phi^{m}_{r} \vec{e}_r +
- N^m (\vec{x}) \phi^{m}_{\theta} \vec{e}_\theta +
- N^m (\vec{x}) \phi^{m}_{z} \vec{e}_z
+ N^m \phi^{m}_{r} \vec{e}_r +
+ N^m \phi^{m}_{\theta} \vec{e}_\theta +
+ N^m \phi^{m}_{z} \vec{e}_z
\bigg)
\Bigg]
\\
&=
\sum\limits_{m=1}^{n_\mathrm{nodes}}
\Bigg[
\phi^{m}_{r } \frac{\partial N^m}{\partial r} \vec{e}_r \vec{e}_r
+ \phi^{m}_{\theta} \frac{\partial N^m}{\partial r} \vec{e}_r \vec{e}_\theta
+ \phi^{m}_{z } \frac{\partial N^m}{\partial r} \vec{e}_r \vec{e}_z
\nonumber
\\
&\quad
- + \phi^{m}_{r } \frac{1}{r} \frac{\partial N^m }{\partial \theta} \vec{e}_\theta \vec{e}_r
- + \phi^{m}_{r } \frac{1}{r} N^m (\vec{x}) \vec{e}_\theta \vec{e}_\theta
- + \phi^{m}_{\theta} \frac{1}{r} \frac{\partial N^m }{\partial \theta} \vec{e}_\theta \vec{e}_\theta
- - \phi^{m}_{\theta} \frac{1}{r} N^m (\vec{x}) \vec{e}_r \vec{e}_\theta
- + \phi^{m}_{z } \frac{1}{r} \frac{\partial N^m }{\partial \theta} \vec{e}_\theta \vec{e}_z
+ + \phi^{m}_{r } \frac{1}{r} \frac{\partial N^m}{\partial \theta} \vec{e}_\theta \vec{e}_r
+ + \phi^{m}_{r } \frac{1}{r} N^m \vec{e}_\theta \vec{e}_\theta
+ + \phi^{m}_{\theta} \frac{1}{r} \frac{\partial N^m}{\partial \theta} \vec{e}_\theta \vec{e}_\theta
+ - \phi^{m}_{\theta} \frac{1}{r} N^m \vec{e}_r \vec{e}_\theta
+ + \phi^{m}_{z } \frac{1}{r} \frac{\partial N^m}{\partial \theta} \vec{e}_\theta \vec{e}_z
\nonumber
\\
&\quad
+ \phi^{m}_{r } \frac{\partial N^m }{\partial z} \vec{e}_z \vec{e}_r
+ \phi^{m}_{\theta} \frac{\partial N^m }{\partial z} \vec{e}_z \vec{e}_\theta
+ \phi^{m}_{z } \frac{\partial N^m }{\partial z} \vec{e}_z \vec{e}_z
\Bigg]
\end{align}
-Which can be denoted in short as
+whereby the spatial dependence was not denoted, but naturally persists. The above can be denoted in short as
\begin{equation}
\vec{\nabla} \vec{\phi} (\vec{x})
\approx
\sum\limits_{m=1}^{n_\mathrm{nodes}}
\mathcal{B}_m (\vec{x}) \cdot \vec{\phi}_m
\end{equation}
where the non-zero components of the third order $\mathcal{B}_m$ read
-\begin{alignat}{3}
- \mathcal{B}^{m}_{r r r } &= \frac{\partial N^m}{\partial r} & \qquad
- \mathcal{B}^{m}_{\theta r r } &= \frac{1}{r} \frac{\partial N^m}{\partial \theta} & \qquad
- \mathcal{B}^{m}_{z r r } &= \frac{\partial N^m}{\partial z}
+\begin{alignat}{4}
+ \mathcal{B}^{m}_{r r r } &= \frac{\partial N^m}{\partial r} & \qquad
+ \mathcal{B}^{m}_{r \theta \theta} &= \frac{\partial N^m}{\partial r} - \frac{1}{r} N^m & \qquad
+ \mathcal{B}^{m}_{r z z } &= \frac{\partial N^m}{\partial r}
\nonumber
\\
- \mathcal{B}^{m}_{r \theta \theta} &= - \frac{1}{r} N^m (\vec{x}) & \qquad
- \mathcal{B}^{m}_{\theta \theta r } &= \frac{1}{r} N^m (\vec{x}) & \qquad
+ \mathcal{B}^{m}_{\theta r r } &= \frac{1}{r} \frac{\partial N^m}{\partial \theta} & \qquad
+ \mathcal{B}^{m}_{\theta \theta \theta} &= \frac{1}{r} \frac{\partial N^m}{\partial \theta} & \qquad
+ \mathcal{B}^{m}_{\theta z z } &= \frac{1}{r} \frac{\partial N^m}{\partial \theta} & \qquad
+ \mathcal{B}^{m}_{\theta \theta r } &= \frac{1}{r} N^m
\nonumber
\\
- \mathcal{B}^{m}_{r \theta \theta} &= \frac{\partial N^m}{\partial r} & \qquad
- \mathcal{B}^{m}_{\theta \theta \theta} &= \frac{1}{r} \frac{\partial N^m}{\partial \theta} & \qquad
- \mathcal{B}^{m}_{z \theta \theta} &= \frac{\partial N^m}{\partial z}
- \\
- \mathcal{B}^{m}_{r z z } &= \frac{\partial N^m}{\partial r} & \qquad
- \mathcal{B}^{m}_{\theta z z } &= \frac{1}{r} \frac{\partial N^m}{\partial \theta} & \qquad
- \mathcal{B}^{m}_{z z z } &= \frac{\partial N^m}{\partial z}
+ \mathcal{B}^{m}_{z r r } &= \frac{\partial N^m}{\partial z} & \qquad
+ \mathcal{B}^{m}_{z \theta \theta} &= \frac{\partial N^m}{\partial z} & \qquad
+ \mathcal{B}^{m}_{z z z } &= \frac{\partial N^m}{\partial z}
\end{alignat}
-Finally a Jacobian $\bm{J}$ can be defined to change coordinates as follows
+\paragraph{Jacobian}
+
+A Jacobian $\bm{J}$ can be defined to change coordinates as follows
\begin{equation}
\vec{\nabla}_\xi = \bm{J} \cdot \vec{\nabla}
\end{equation}
which in components reads
\begin{equation}
\vec{e}_\xi \frac{\partial}{\partial \xi } +
\vec{e}_\eta \frac{\partial}{\partial \eta} +
\vec{e}_\mu \frac{\partial}{\partial \mu}
=
\bm{J} \cdot
\left(
\vec{e}_r \frac{\partial}{\partial r } +
\frac{1}{r} \vec{e}_\theta \frac{\partial}{\partial \theta} +
\vec{e}_z \frac{\partial}{\partial z }
\right)
\end{equation}
The Jacobian is then obtained as
\begin{align}
\bm{J}
&=
\vec{e}_\xi \frac{\partial r }{\partial \xi } \vec{e}_r +
\vec{e}_\xi r \frac{\partial \theta}{\partial \xi } \vec{e}_\theta +
\vec{e}_\xi \frac{\partial z }{\partial \xi } \vec{e}_z
\nonumber
\\
&+
\vec{e}_\eta \frac{\partial r }{\partial \eta} \vec{e}_r +
\vec{e}_\eta r \frac{\partial \theta}{\partial \eta} \vec{e}_\theta +
\vec{e}_\eta \frac{\partial z }{\partial \eta} \vec{e}_z
\nonumber
\\
&+
\vec{e}_\mu \frac{\partial r }{\partial \mu } \vec{e}_r +
\vec{e}_\mu r \frac{\partial \theta}{\partial \mu } \vec{e}_\theta +
\vec{e}_\mu \frac{\partial z }{\partial \mu } \vec{e}_z
\end{align}
+\paragraph{Axisymmetric}
+
+Assuming axisymmetry implies assuming that there can be no variations along the $\vec{e}_\theta$ direction. This implies that any derivate $\partial / \partial \theta = 0$ and consequently that there are no degrees-of-freedom in the $\theta$ direction, i.e.\ $u_\theta = 0$. This implies that the only relevant non-zero components of the third order $\mathcal{B}_m$ read
+\begin{alignat}{4}
+ \mathcal{B}^{m}_{r r r } &= \frac{\partial N^m}{\partial r} & \qquad
+ \mathcal{B}^{m}_{r z z } &= \frac{\partial N^m}{\partial r}
+ \nonumber
+ \\
+ \mathcal{B}^{m}_{\theta \theta r } &= \frac{1}{r} N^m
+ \nonumber
+ \\
+ \mathcal{B}^{m}_{z r r } &= \frac{\partial N^m}{\partial z} & \qquad
+ \mathcal{B}^{m}_{z z z } &= \frac{\partial N^m}{\partial z}
+\end{alignat}
+The Jacobian is furthermore reduced to
+\begin{align}
+ \bm{J}
+ &=
+ \vec{e}_\xi \frac{\partial r }{\partial \xi } \vec{e}_r +
+ \vec{e}_\xi \frac{\partial z }{\partial \xi } \vec{e}_z
+ \nonumber
+ \\
+ &+
+ \vec{e}_\eta \frac{\partial r }{\partial \eta} \vec{e}_r +
+ \vec{e}_\eta \frac{\partial z }{\partial \eta} \vec{e}_z
+\end{align}
% \bibliography{example_refs}
\end{document}
diff --git a/docs/theory/symbolic/axisymmetric.py b/docs/theory/symbolic/axisymmetric.py
new file mode 100644
index 0000000..f5f6a3b
--- /dev/null
+++ b/docs/theory/symbolic/axisymmetric.py
@@ -0,0 +1,42 @@
+
+import numpy as np
+
+B = np.zeros((3,3,3))
+BT = np.zeros((3,3,3))
+
+B[0,1,1] = 1. # B(r \theta \theta )
+B[1,1,0] = 1. # B(\theta \theta r )
+B[0,0,0] = 1. # B(r r r )
+B[0,1,1] = 1. # B(r \theta \theta )
+B[0,2,2] = 1. # B(r z z )
+B[2,0,0] = 1. # B(z r r )
+B[2,1,1] = 1. # B(z \theta \theta )
+B[2,2,2] = 1. # B(z z z )
+
+for i in range(3):
+ for j in range(3):
+ for k in range(3):
+ BT[k,j,i] = B[i,j,k]
+
+C = np.ones((3,3,3,3))
+
+K = [['' for j in range(3)] for i in range(3)]
+
+def st(i):
+ if i == 0: return 'r'
+ if i == 1: return 't'
+ if i == 2: return 'z'
+
+for i in range(3):
+ for j in range(3):
+ for k in range(3):
+ for l in range(3):
+ for m in range(3):
+ for n in range(3):
+ if BT[i,j,k] * C[k,j,l,m] * B[m,l,n]:
+ K[i][n] += ' + ' + 'B(%s,%s,%s)'%(st(k),st(j),st(i)) + '*C(%s,%s,%s,%s)*'%(st(k),st(j),st(l),st(m)) + 'B(%s,%s,%s)'%(st(m),st(l),st(n))
+
+i = 0; j = 0; print(i,j,K[i][j])
+i = 0; j = 2; print(i,j,K[i][j])
+i = 2; j = 0; print(i,j,K[i][j])
+i = 2; j = 2; print(i,j,K[i][j])
diff --git a/docs/theory/unsrtnat.bst b/docs/theory/unsrtnat.bst
index 1c7eddf..a212b66 100644
--- a/docs/theory/unsrtnat.bst
+++ b/docs/theory/unsrtnat.bst
@@ -1,1343 +1,1343 @@
%% File: `unsrtnat.bst'
-%% A modification of `unsrt.bst' for use with natbib package
+%% A modification of `unsrt.bst' for use with natbib package
%%
%% Copyright 1993-2007 Patrick W Daly
%% Max-Planck-Institut f\"ur Sonnensystemforschung
%% Max-Planck-Str. 2
%% D-37191 Katlenburg-Lindau
%% Germany
%% E-mail: daly@mps.mpg.de
%%
%% This program can be redistributed and/or modified under the terms
%% of the LaTeX Project Public License Distributed from CTAN
%% archives in directory macros/latex/base/lppl.txt; either
%% version 1 of the License, or any later version.
%%
% Version and source file information:
% \ProvidesFile{natbst.mbs}[2007/11/26 1.93 (PWD)]
%
% BibTeX `plainnat' family
% version 0.99b for BibTeX versions 0.99a or later,
% for LaTeX versions 2.09 and 2e.
%
% For use with the `natbib.sty' package; emulates the corresponding
% member of the `plain' family, but with author-year citations.
%
% With version 6.0 of `natbib.sty', it may also be used for numerical
% citations, while retaining the commands \citeauthor, \citefullauthor,
% and \citeyear to print the corresponding information.
%
% For version 7.0 of `natbib.sty', the KEY field replaces missing
% authors/editors, and the date is left blank in \bibitem.
%
% Includes field EID for the sequence/citation number of electronic journals
% which is used instead of page numbers.
%
% Includes fields ISBN and ISSN.
%
% Includes field URL for Internet addresses.
%
% Includes field DOI for Digital Object Idenfifiers.
%
% Works best with the url.sty package of Donald Arseneau.
%
% Works with identical authors and year are further sorted by
% title text, as in the standard plain.bst etc.
%
ENTRY
{ address
author
booktitle
chapter
doi
- eprint
+ arxivid
eid
edition
editor
howpublished
institution
isbn
issn
journal
key
month
note
number
organization
pages
publisher
school
series
title
type
url
volume
year
}
{}
{ label extra.label sort.label short.list }
INTEGERS { output.state before.all mid.sentence after.sentence after.block }
FUNCTION {init.state.consts}
{ #0 'before.all :=
#1 'mid.sentence :=
#2 'after.sentence :=
#3 'after.block :=
}
STRINGS { s t }
FUNCTION {output.nonnull}
{ 's :=
output.state mid.sentence =
{ ", " * write$ }
{ output.state after.block =
{ add.period$ write$
newline$
"\newblock " write$
}
{ output.state before.all =
'write$
{ add.period$ " " * write$ }
if$
}
if$
mid.sentence 'output.state :=
}
if$
s
}
FUNCTION {output}
{ duplicate$ empty$
'pop$
'output.nonnull
if$
}
FUNCTION {output.check}
{ 't :=
duplicate$ empty$
{ pop$ "empty " t * " in " * cite$ * warning$ }
'output.nonnull
if$
}
FUNCTION {fin.entry}
{ add.period$
write$
newline$
}
FUNCTION {new.block}
{ output.state before.all =
'skip$
{ after.block 'output.state := }
if$
}
FUNCTION {new.sentence}
{ output.state after.block =
'skip$
{ output.state before.all =
'skip$
{ after.sentence 'output.state := }
if$
}
if$
}
FUNCTION {not}
{ { #0 }
{ #1 }
if$
}
FUNCTION {and}
{ 'skip$
{ pop$ #0 }
if$
}
FUNCTION {or}
{ { pop$ #1 }
'skip$
if$
}
FUNCTION {new.block.checka}
{ empty$
'skip$
'new.block
if$
}
FUNCTION {new.block.checkb}
{ empty$
swap$ empty$
and
'skip$
'new.block
if$
}
FUNCTION {new.sentence.checka}
{ empty$
'skip$
'new.sentence
if$
}
FUNCTION {new.sentence.checkb}
{ empty$
swap$ empty$
and
'skip$
'new.sentence
if$
}
FUNCTION {field.or.null}
{ duplicate$ empty$
{ pop$ "" }
'skip$
if$
}
FUNCTION {emphasize}
{ duplicate$ empty$
{ pop$ "" }
{ "\emph{" swap$ * "}" * }
if$
}
INTEGERS { nameptr namesleft numnames }
FUNCTION {format.names}
{ 's :=
#1 'nameptr :=
s num.names$ 'numnames :=
numnames 'namesleft :=
{ namesleft #0 > }
{ s nameptr "{ff~}{vv~}{ll}{, jj}" format.name$ 't :=
nameptr #1 >
{ namesleft #1 >
{ ", " * t * }
{ numnames #2 >
{ "," * }
'skip$
if$
t "others" =
{ " et~al." * }
{ " and " * t * }
if$
}
if$
}
't
if$
nameptr #1 + 'nameptr :=
namesleft #1 - 'namesleft :=
}
while$
}
FUNCTION {format.key}
{ empty$
{ key field.or.null }
{ "" }
if$
}
FUNCTION {format.authors}
{ author empty$
{ "" }
{ author format.names }
if$
}
FUNCTION {format.editors}
{ editor empty$
{ "" }
{ editor format.names
editor num.names$ #1 >
{ ", editors" * }
{ ", editor" * }
if$
}
if$
}
FUNCTION {format.isbn}
{ isbn empty$
{ "" }
{ new.block "ISBN " isbn * }
if$
}
FUNCTION {format.issn}
{ issn empty$
{ "" }
{ new.block "ISSN " issn * }
if$
}
FUNCTION {format.url}
{ url empty$
{ "" }
{ new.block "URL \url{" url * "}" * }
if$
}
FUNCTION {format.doi}
{ doi empty$
{ "" }
{ new.block "\doi{" doi * "}" * }
if$
}
-FUNCTION {format.eprint}
-{ eprint empty$
+FUNCTION {format.arxivid}
+{ arxivid empty$
{ "" }
- { new.block "\eprint{" eprint * "}" * }
+ { new.block "\arxivid{" arxivid * "}" * }
if$
}
FUNCTION {format.title}
{ title empty$
{ "" }
{ title "t" change.case$ }
if$
}
FUNCTION {format.full.names}
{'s :=
#1 'nameptr :=
s num.names$ 'numnames :=
numnames 'namesleft :=
{ namesleft #0 > }
{ s nameptr
"{vv~}{ll}" format.name$ 't :=
nameptr #1 >
{
namesleft #1 >
{ ", " * t * }
{
numnames #2 >
{ "," * }
'skip$
if$
t "others" =
{ " et~al." * }
{ " and " * t * }
if$
}
if$
}
't
if$
nameptr #1 + 'nameptr :=
namesleft #1 - 'namesleft :=
}
while$
}
FUNCTION {author.editor.full}
{ author empty$
{ editor empty$
{ "" }
{ editor format.full.names }
if$
}
{ author format.full.names }
if$
}
FUNCTION {author.full}
{ author empty$
{ "" }
{ author format.full.names }
if$
}
FUNCTION {editor.full}
{ editor empty$
{ "" }
{ editor format.full.names }
if$
}
FUNCTION {make.full.names}
{ type$ "book" =
type$ "inbook" =
or
'author.editor.full
{ type$ "proceedings" =
'editor.full
'author.full
if$
}
if$
}
FUNCTION {output.bibitem}
{ newline$
"\bibitem[" write$
label write$
")" make.full.names duplicate$ short.list =
{ pop$ }
{ * }
if$
"]{" * write$
cite$ write$
"}" write$
newline$
""
before.all 'output.state :=
}
FUNCTION {n.dashify}
{ 't :=
""
{ t empty$ not }
{ t #1 #1 substring$ "-" =
{ t #1 #2 substring$ "--" = not
{ "--" *
t #2 global.max$ substring$ 't :=
}
{ { t #1 #1 substring$ "-" = }
{ "-" *
t #2 global.max$ substring$ 't :=
}
while$
}
if$
}
{ t #1 #1 substring$ *
t #2 global.max$ substring$ 't :=
}
if$
}
while$
}
FUNCTION {format.date}
{ year duplicate$ empty$
{ "empty year in " cite$ * warning$
pop$ "" }
'skip$
if$
month empty$
'skip$
{ month
" " * swap$ *
}
if$
extra.label *
}
FUNCTION {format.btitle}
{ title emphasize
}
FUNCTION {tie.or.space.connect}
{ duplicate$ text.length$ #3 <
{ "~" }
{ " " }
if$
swap$ * *
}
FUNCTION {either.or.check}
{ empty$
'pop$
{ "can't use both " swap$ * " fields in " * cite$ * warning$ }
if$
}
FUNCTION {format.bvolume}
{ volume empty$
{ "" }
{ "volume" volume tie.or.space.connect
series empty$
'skip$
{ " of " * series emphasize * }
if$
"volume and number" number either.or.check
}
if$
}
FUNCTION {format.number.series}
{ volume empty$
{ number empty$
{ series field.or.null }
{ output.state mid.sentence =
{ "number" }
{ "Number" }
if$
number tie.or.space.connect
series empty$
{ "there's a number but no series in " cite$ * warning$ }
{ " in " * series * }
if$
}
if$
}
{ "" }
if$
}
FUNCTION {format.edition}
{ edition empty$
{ "" }
{ output.state mid.sentence =
{ edition "l" change.case$ " edition" * }
{ edition "t" change.case$ " edition" * }
if$
}
if$
}
INTEGERS { multiresult }
FUNCTION {multi.page.check}
{ 't :=
#0 'multiresult :=
{ multiresult not
t empty$ not
and
}
{ t #1 #1 substring$
duplicate$ "-" =
swap$ duplicate$ "," =
swap$ "+" =
or or
{ #1 'multiresult := }
{ t #2 global.max$ substring$ 't := }
if$
}
while$
multiresult
}
FUNCTION {format.pages}
{ pages empty$
{ "" }
{ pages multi.page.check
{ "pages" pages n.dashify tie.or.space.connect }
{ "page" pages tie.or.space.connect }
if$
}
if$
}
FUNCTION {format.eid}
{ eid empty$
{ "" }
{ "art." eid tie.or.space.connect }
if$
}
FUNCTION {format.vol.num.pages}
{ volume field.or.null
number empty$
'skip$
{ "\penalty0 (" number * ")" * *
volume empty$
{ "there's a number but no volume in " cite$ * warning$ }
'skip$
if$
}
if$
pages empty$
'skip$
{ duplicate$ empty$
{ pop$ format.pages }
{ ":\penalty0 " * pages n.dashify * }
if$
}
if$
}
FUNCTION {format.vol.num.eid}
{ volume field.or.null
number empty$
'skip$
{ "\penalty0 (" number * ")" * *
volume empty$
{ "there's a number but no volume in " cite$ * warning$ }
'skip$
if$
}
if$
eid empty$
'skip$
{ duplicate$ empty$
{ pop$ format.eid }
{ ":\penalty0 " * eid * }
if$
}
if$
}
FUNCTION {format.chapter.pages}
{ chapter empty$
'format.pages
{ type empty$
{ "chapter" }
{ type "l" change.case$ }
if$
chapter tie.or.space.connect
pages empty$
'skip$
{ ", " * format.pages * }
if$
}
if$
}
FUNCTION {format.in.ed.booktitle}
{ booktitle empty$
{ "" }
{ editor empty$
{ "In " booktitle emphasize * }
{ "In " format.editors * ", " * booktitle emphasize * }
if$
}
if$
}
FUNCTION {empty.misc.check}
{ author empty$ title empty$ howpublished empty$
month empty$ year empty$ note empty$
and and and and and
key empty$ not and
{ "all relevant fields are empty in " cite$ * warning$ }
'skip$
if$
}
FUNCTION {format.thesis.type}
{ type empty$
'skip$
{ pop$
type "t" change.case$
}
if$
}
FUNCTION {format.tr.number}
{ type empty$
{ "Technical Report" }
'type
if$
number empty$
{ "t" change.case$ }
{ number tie.or.space.connect }
if$
}
FUNCTION {format.article.crossref}
{ key empty$
{ journal empty$
{ "need key or journal for " cite$ * " to crossref " * crossref *
warning$
""
}
{ "In \emph{" journal * "}" * }
if$
}
{ "In " }
if$
" \citet{" * crossref * "}" *
}
FUNCTION {format.book.crossref}
{ volume empty$
{ "empty volume in " cite$ * "'s crossref of " * crossref * warning$
"In "
}
{ "Volume" volume tie.or.space.connect
" of " *
}
if$
editor empty$
editor field.or.null author field.or.null =
or
{ key empty$
{ series empty$
{ "need editor, key, or series for " cite$ * " to crossref " *
crossref * warning$
"" *
}
{ "\emph{" * series * "}" * }
if$
}
'skip$
if$
}
'skip$
if$
" \citet{" * crossref * "}" *
}
FUNCTION {format.incoll.inproc.crossref}
{ editor empty$
editor field.or.null author field.or.null =
or
{ key empty$
{ booktitle empty$
{ "need editor, key, or booktitle for " cite$ * " to crossref " *
crossref * warning$
""
}
{ "In \emph{" booktitle * "}" * }
if$
}
{ "In " }
if$
}
{ "In " }
if$
" \citet{" * crossref * "}" *
}
FUNCTION {article}
{ output.bibitem
format.authors "author" output.check
author format.key output
new.block
format.title "title" output.check
new.block
crossref missing$
{ journal emphasize "journal" output.check
eid empty$
{ format.vol.num.pages output }
{ format.vol.num.eid output }
if$
format.date "year" output.check
}
{ format.article.crossref output.nonnull
eid empty$
{ format.pages output }
{ format.eid output }
if$
}
if$
format.issn output
format.doi output
- format.eprint output
+ format.arxivid output
format.url output
new.block
note output
fin.entry
}
FUNCTION {book}
{ output.bibitem
author empty$
{ format.editors "author and editor" output.check
editor format.key output
}
{ format.authors output.nonnull
crossref missing$
{ "author and editor" editor either.or.check }
'skip$
if$
}
if$
new.block
format.btitle "title" output.check
crossref missing$
{ format.bvolume output
new.block
format.number.series output
new.sentence
publisher "publisher" output.check
address output
}
{ new.block
format.book.crossref output.nonnull
}
if$
format.edition output
format.date "year" output.check
format.isbn output
format.doi output
- format.eprint output
+ format.arxivid output
format.url output
new.block
note output
fin.entry
}
FUNCTION {booklet}
{ output.bibitem
format.authors output
author format.key output
new.block
format.title "title" output.check
howpublished address new.block.checkb
howpublished output
address output
format.date output
format.isbn output
format.doi output
- format.eprint output
+ format.arxivid output
format.url output
new.block
note output
fin.entry
}
FUNCTION {inbook}
{ output.bibitem
author empty$
{ format.editors "author and editor" output.check
editor format.key output
}
{ format.authors output.nonnull
crossref missing$
{ "author and editor" editor either.or.check }
'skip$
if$
}
if$
new.block
format.btitle "title" output.check
crossref missing$
{ format.bvolume output
format.chapter.pages "chapter and pages" output.check
new.block
format.number.series output
new.sentence
publisher "publisher" output.check
address output
}
{ format.chapter.pages "chapter and pages" output.check
new.block
format.book.crossref output.nonnull
}
if$
format.edition output
format.date "year" output.check
format.isbn output
format.doi output
- format.eprint output
+ format.arxivid output
format.url output
new.block
note output
fin.entry
}
FUNCTION {incollection}
{ output.bibitem
format.authors "author" output.check
author format.key output
new.block
format.title "title" output.check
new.block
crossref missing$
{ format.in.ed.booktitle "booktitle" output.check
format.bvolume output
format.number.series output
format.chapter.pages output
new.sentence
publisher "publisher" output.check
address output
format.edition output
format.date "year" output.check
}
{ format.incoll.inproc.crossref output.nonnull
format.chapter.pages output
}
if$
format.isbn output
format.doi output
- format.eprint output
+ format.arxivid output
format.url output
new.block
note output
fin.entry
}
FUNCTION {inproceedings}
{ output.bibitem
format.authors "author" output.check
author format.key output
new.block
format.title "title" output.check
new.block
crossref missing$
{ format.in.ed.booktitle "booktitle" output.check
format.bvolume output
format.number.series output
format.pages output
address empty$
{ organization publisher new.sentence.checkb
organization output
publisher output
format.date "year" output.check
}
{ address output.nonnull
format.date "year" output.check
new.sentence
organization output
publisher output
}
if$
}
{ format.incoll.inproc.crossref output.nonnull
format.pages output
}
if$
format.isbn output
format.doi output
- format.eprint output
+ format.arxivid output
format.url output
new.block
note output
fin.entry
}
FUNCTION {conference} { inproceedings }
FUNCTION {manual}
{ output.bibitem
format.authors output
author format.key output
new.block
format.btitle "title" output.check
organization address new.block.checkb
organization output
address output
format.edition output
format.date output
format.url output
new.block
note output
fin.entry
}
FUNCTION {mastersthesis}
{ output.bibitem
format.authors "author" output.check
author format.key output
new.block
format.title "title" output.check
new.block
"Master's thesis" format.thesis.type output.nonnull
school "school" output.check
address output
format.date "year" output.check
format.url output
new.block
note output
fin.entry
}
FUNCTION {misc}
{ output.bibitem
format.authors output
author format.key output
title howpublished new.block.checkb
format.title output
howpublished new.block.checka
howpublished output
format.date output
format.issn output
format.url output
new.block
note output
fin.entry
empty.misc.check
}
FUNCTION {phdthesis}
{ output.bibitem
format.authors "author" output.check
author format.key output
new.block
format.btitle "title" output.check
new.block
"PhD thesis" format.thesis.type output.nonnull
school "school" output.check
address output
format.date "year" output.check
format.url output
new.block
note output
fin.entry
}
FUNCTION {proceedings}
{ output.bibitem
format.editors output
editor format.key output
new.block
format.btitle "title" output.check
format.bvolume output
format.number.series output
address output
format.date "year" output.check
new.sentence
organization output
publisher output
format.isbn output
format.doi output
- format.eprint output
+ format.arxivid output
format.url output
new.block
note output
fin.entry
}
FUNCTION {techreport}
{ output.bibitem
format.authors "author" output.check
author format.key output
new.block
format.title "title" output.check
new.block
format.tr.number output.nonnull
institution "institution" output.check
address output
format.date "year" output.check
format.url output
new.block
note output
fin.entry
}
FUNCTION {unpublished}
{ output.bibitem
format.authors "author" output.check
author format.key output
new.block
format.title "title" output.check
new.block
note "note" output.check
format.date output
format.url output
fin.entry
}
FUNCTION {default.type} { misc }
MACRO {jan} {"January"}
MACRO {feb} {"February"}
MACRO {mar} {"March"}
MACRO {apr} {"April"}
MACRO {may} {"May"}
MACRO {jun} {"June"}
MACRO {jul} {"July"}
MACRO {aug} {"August"}
MACRO {sep} {"September"}
MACRO {oct} {"October"}
MACRO {nov} {"November"}
MACRO {dec} {"December"}
MACRO {acmcs} {"ACM Computing Surveys"}
MACRO {acta} {"Acta Informatica"}
MACRO {cacm} {"Communications of the ACM"}
MACRO {ibmjrd} {"IBM Journal of Research and Development"}
MACRO {ibmsj} {"IBM Systems Journal"}
MACRO {ieeese} {"IEEE Transactions on Software Engineering"}
MACRO {ieeetc} {"IEEE Transactions on Computers"}
MACRO {ieeetcad}
{"IEEE Transactions on Computer-Aided Design of Integrated Circuits"}
MACRO {ipl} {"Information Processing Letters"}
MACRO {jacm} {"Journal of the ACM"}
MACRO {jcss} {"Journal of Computer and System Sciences"}
MACRO {scp} {"Science of Computer Programming"}
MACRO {sicomp} {"SIAM Journal on Computing"}
MACRO {tocs} {"ACM Transactions on Computer Systems"}
MACRO {tods} {"ACM Transactions on Database Systems"}
MACRO {tog} {"ACM Transactions on Graphics"}
MACRO {toms} {"ACM Transactions on Mathematical Software"}
MACRO {toois} {"ACM Transactions on Office Information Systems"}
MACRO {toplas} {"ACM Transactions on Programming Languages and Systems"}
MACRO {tcs} {"Theoretical Computer Science"}
READ
FUNCTION {sortify}
{ purify$
"l" change.case$
}
INTEGERS { len }
FUNCTION {chop.word}
{ 's :=
'len :=
s #1 len substring$ =
{ s len #1 + global.max$ substring$ }
's
if$
}
FUNCTION {format.lab.names}
{ 's :=
s #1 "{vv~}{ll}" format.name$
s num.names$ duplicate$
#2 >
{ pop$ " et~al." * }
{ #2 <
'skip$
{ s #2 "{ff }{vv }{ll}{ jj}" format.name$ "others" =
{ " et~al." * }
{ " and " * s #2 "{vv~}{ll}" format.name$ * }
if$
}
if$
}
if$
}
FUNCTION {author.key.label}
{ author empty$
{ key empty$
{ cite$ #1 #3 substring$ }
'key
if$
}
{ author format.lab.names }
if$
}
FUNCTION {author.editor.key.label}
{ author empty$
{ editor empty$
{ key empty$
{ cite$ #1 #3 substring$ }
'key
if$
}
{ editor format.lab.names }
if$
}
{ author format.lab.names }
if$
}
FUNCTION {author.key.organization.label}
{ author empty$
{ key empty$
{ organization empty$
{ cite$ #1 #3 substring$ }
{ "The " #4 organization chop.word #3 text.prefix$ }
if$
}
'key
if$
}
{ author format.lab.names }
if$
}
FUNCTION {editor.key.organization.label}
{ editor empty$
{ key empty$
{ organization empty$
{ cite$ #1 #3 substring$ }
{ "The " #4 organization chop.word #3 text.prefix$ }
if$
}
'key
if$
}
{ editor format.lab.names }
if$
}
FUNCTION {calc.short.authors}
{ type$ "book" =
type$ "inbook" =
or
'author.editor.key.label
{ type$ "proceedings" =
'editor.key.organization.label
{ type$ "manual" =
'author.key.organization.label
'author.key.label
if$
}
if$
}
if$
'short.list :=
}
FUNCTION {calc.label}
{ calc.short.authors
short.list
"("
*
year duplicate$ empty$
short.list key field.or.null = or
{ pop$ "" }
'skip$
if$
*
'label :=
}
INTEGERS { seq.num }
FUNCTION {init.seq}
{ #0 'seq.num :=}
EXECUTE {init.seq}
FUNCTION {int.to.fix}
{ "000000000" swap$ int.to.str$ *
#-1 #10 substring$
}
FUNCTION {presort}
{ calc.label
label sortify
" "
*
seq.num #1 + 'seq.num :=
seq.num int.to.fix
'sort.label :=
sort.label *
#1 entry.max$ substring$
'sort.key$ :=
}
ITERATE {presort}
SORT
STRINGS { longest.label last.label next.extra }
INTEGERS { longest.label.width last.extra.num number.label }
FUNCTION {initialize.longest.label}
{ "" 'longest.label :=
#0 int.to.chr$ 'last.label :=
"" 'next.extra :=
#0 'longest.label.width :=
#0 'last.extra.num :=
#0 'number.label :=
}
FUNCTION {forward.pass}
{ last.label label =
{ last.extra.num #1 + 'last.extra.num :=
last.extra.num int.to.chr$ 'extra.label :=
}
{ "a" chr.to.int$ 'last.extra.num :=
"" 'extra.label :=
label 'last.label :=
}
if$
number.label #1 + 'number.label :=
}
FUNCTION {reverse.pass}
{ next.extra "b" =
{ "a" 'extra.label := }
'skip$
if$
extra.label 'next.extra :=
extra.label
duplicate$ empty$
'skip$
{ "{\natexlab{" swap$ * "}}" * }
if$
'extra.label :=
label extra.label * 'label :=
}
EXECUTE {initialize.longest.label}
ITERATE {forward.pass}
REVERSE {reverse.pass}
FUNCTION {bib.sort.order}
{ sort.label 'sort.key$ :=
}
ITERATE {bib.sort.order}
SORT
FUNCTION {begin.bib}
{ preamble$ empty$
'skip$
{ preamble$ write$ newline$ }
if$
"\begin{thebibliography}{" number.label int.to.str$ * "}" *
write$ newline$
"\providecommand{\natexlab}[1]{#1}"
write$ newline$
"\providecommand{\url}[1]{\texttt{#1}}"
write$ newline$
"\expandafter\ifx\csname urlstyle\endcsname\relax"
write$ newline$
- " \providecommand{\eprint}[1]{eprint: #1}\else"
+ " \providecommand{\arxivid}[1]{arxivid: #1}\else"
write$ newline$
- " \providecommand{\eprint}{eprint: \begingroup \urlstyle{rm}\Url}\fi"
+ " \providecommand{\arxivid}{arxivid: \begingroup \urlstyle{rm}\Url}\fi"
write$ newline$
"\expandafter\ifx\csname urlstyle\endcsname\relax"
write$ newline$
" \providecommand{\doi}[1]{doi: #1}\else"
write$ newline$
" \providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi"
write$ newline$
}
EXECUTE {begin.bib}
EXECUTE {init.state.consts}
ITERATE {call.type$}
FUNCTION {end.bib}
{ newline$
"\end{thebibliography}" write$ newline$
}
EXECUTE {end.bib}